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ABSTRACT 


Vincent (2006) developed a technique for remotely measuring Aerosol Optical Depth 
(AOD) using commercial high-resolution satellite imagery. This technique measured the 
radiance difference between a building shadow and an adjacent sunlit region with the 
same surface reflectance to calculate Total Optical Depth (TOD). AOD is then 
determined by subtracting Rayleigh scattering from TOD. The procedure for making this 
calculation was time consuming, particularly locating suitable shadows within the region 
of interest. This paper outlines a fully automated method of performing the AOD 
calculation and examining shadow properties. The automated method relies on high- 
resolution Digital Surface Model (DSM) data collected using a Light Detection and 
Ranging (LIDAR) sensor coupled with sun and satellite geometry to identify shadow 
regions. Configuration settings allowed for specific regions in the shadow and sunlit area 
to be selected before determining their respective radiances. Finally, a technique for 
aligning the satellite and DSM pixels was developed to correct for small differences 
between the datasets. Results from the automated method were compared with 
AERONET data for validation. The automated method using WorldView-1 and 
QuickBird imagery worked best at Solar Village, Saudi Arabia and an area northeast of 
Washington, D.C., which included the Goddard Space Flight Center. Testing of 
IKONOS multispectral imagery suggested the resolution is inadequate in urban settings. 
Testing in areas that included downtown regions in Houston, TX and Baltimore, MD 
identified weaknesses in the alignment algorithm. 
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TOD 

Total Optical Depth 
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I. 


INTRODUCTION 


Aerosols play an important part in the earth’s radiation budget directly (through 
scattering and absorption of incoming solar radiation and outgoing longwave radiation) 
and indirectly (by affecting cloud properties) (IPCC 2007). The overall impact depends 
on the aerosol optical properties that vary with wavelength, relative humidity, and the 
horizontal and vertical distribution, all of these varying with time (Hayward and Boucher 
2000, IPCC 2001). Although aerosols can contribute positively to global heating, their 
overall effect is a net cooling effect. This cooling counteracts 25-50% of the heating 
from anthropogenic greenhouse gases (IPCC 2007). Similarly, aerosols significantly 
attenuate transmission of electromagnetic radiation in the visible and near-infrared (NIR) 
through scattering and absorption, limiting the effectiveness of intelligence collection and 
weapon systems that operate in those portions of the spectrum by reducing the amount of 
radiation reaching the sensors. Aerosol optical depth (AOD) is a measure of the opacity 
of the atmosphere due to aerosols. It is wavelength-dependent and is directly 
proportional to the size and concentration of the aerosol particles present. Accurately 
characterizing the radiative forcing due to aerosols is important when running climate 
models (IPCC 2007). However, the significant spatial and temporal variability makes it 
difficult to directly measure aerosol properties (size and concentration) in the atmosphere 
globally. As a result, numerous algorithms using satellite imagery to determine aerosol 
properties have been developed. 

Most satellite techniques for characterizing aerosols with visible and NIR 
wavelength sensors rely on surfaces with low reflectance in the visible channels or over 
surfaces where detailed knowledge of the surface reflectance exists. A few satellites 
(Polarization and Directionality of the Earth’s Reflectance - POLDER and Multi-angle 
Imaging Spectro-Radiometer - MISR) have dedicated sensors for measuring aerosol 
properties (IPCC 2007). These sensors use multiple look angles to model the atmosphere 
but are limited to low surface reflectance regions (Veefkind et al. 1998, Martonchik et al. 
2004, Deuze et al. 2000). Vincent (2006) developed a technique to characterize aerosols 
using commercial high-resolution satellite data in the visible and NIR portions of the 
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spectrum, over any terrain, without prior knowledge of or constraints on surface 
reflectance. This technique, henceforth referred to as the shadow technique or method, 
compares the radiance inside the shadow with radiance just outside the shadow (sunlit) to 
calculate AOD. Although promising, Vincent’s method was labor intensive, making 
large case studies and practical application difficult. Since that time, additional research 
has refined the technique and sped up the process (Rivenbark 2009, Sweat 2008, 
Dombrock 2007, Evans 2007). However, those efforts stopped well short of full 
automation. 

This research develops and tests a method called the Fully Automated Shadow 
Technique, or FASTEC, of fully automating the shadow technique for any sensor or 
wavelength provided a high resolution surface elevation map is available. In addition to 
identifying and calculating AOD for specific shadows or every usable shadow in an 
image, FASTEC enables detailed examination of shadow and sunlit region characteristics 
for determining pixel selection criteria. Chapter II describes the shadow technique as 
developed by Vincent (2006). It also lists and briefly describes other methods of 
measuring AOD. Chapter III outlines the automated approach and its implementation. 
Chapter IV details the results of testing the automated method. Chapter V discusses 
conclusions and offers suggestions for future research with the automated method. 
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II. BACKGROUND 


A. OTHER SATELLITE-BASED AOD RETRIEVAL METHODS 

Numerous satellite-based techniques for measuring AOD exist. These techniques 
typically work with specific wavelength bands, minimum horizontal resolutions, and 
surface reflectance ranges. Figure 1 is a global seasonal average composite of AOD from 
the Moderate Resolution Imaging Spectrometer (MODIS) sensor, demonstrating the high 
spatial and seasonal variability of AOD globally. Additionally, Figure 1 shows large 
areas where AOD is not calculated because existing MODIS techniques do not work over 
desert or snow-covered regions due to high surface reflectance. 


January to March 2001 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 

Total Aerosol Optical Depth 

Figure 1. Aerosol optical depth at 0.55 /um using the MODIS sensor (from IPCC 2007). 
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Surface reflectance is negligible over water for near-infrared (NIR) wavelengths 
and AOD is calculated using the dark object method, which is described in Chapter II.A.2 
(von Hoyningen-Huene 2003). Here total optical depth (TOD) is computed by treating 
the low surface reflectance as negligible and assuming all radiance received by the 
satellite is scattered spaceward by aerosol particles and gases. After removing the 
contribution from Rayleigh scattering (depends on surface pressure), the contribution 
from aerosol scattering dominates the remaining signal. Over land, determining optical 
depth is less straightforward since the contribution from surface reflectance is non-zero 
and varies significantly with surface type and wavelength. To calculate AOD over land, 
more complex methods like contrast reduction, dark object, and multi-angle were 
developed, and these are described as follows: 

1. Contrast Reduction Method 

The contrast reduction method computes optical depth from the measured 
radiance difference between a distinct surface and its surrounding area (Odell and 
Weinman 1975). It requires prior knowledge of the two surface albedos, sun and 
satellite zenith angles, and lookup tables based on the mean surface albedo, the viewing 
zenith angle, and a contrast transmission function to compute optical depth. As defined 
by Odell and Weinman (1975), the contrast transmission function is proportional to the 
difference between the reflected radiance off the surfaces and inversely proportional to 
the difference between their surface albedos. Odell and Weinman (1975) showed that it 
is possible to exploit shadows cast by large objects to detennine optical depth without 
prior knowledge of the surface features. This forms the basis of the shadow method 
developed by Vincent (2006). Later, Kaufman and Joseph (1982) applied a similar 
approach for regions where the surface reflectance undergoes a step-wise change, such as 
along a shoreline or forest border, to compute optical depth and surface reflectivity. This 
eliminates the previous requirement that surface reflectance be known. Additionally, 
using two different surface regions allows this approach to be independent of satellite 
calibration and produces more accurate results than earlier contrast reduction methods. 
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Dark Object Method 


The dark object method is usable for wavelengths with low surface reflectance. 
In such locations it is possible to estimate the small portion of the top-of-atmosphere 
(TOA) signal, which is due to surface reflection. Kaufman and Sendra (1988) detail a 
rapid, automated algorithm to compute optical thickness, using the dark object method 
over regions covered by dark vegetation. Their approach assumes the surface reflectance 
of those regions to estimate optical thickness and was accurate to within 0.20 for band 1 
(0.5-0.6 pm) and 0.05 for band 2 (0.6-0.7 pm). The dark object method is ineffective 
over surfaces with high surface reflectance, such as snow or sand. Other dark object 
algorithms, such as that described by von Hoyningen-Huene et al. (2003), employ the 
same fundamental idea. 

3. Multi-angle method 

Multi-angle methods use special sensors that examine a region from multiple look 
angles. Sensors such as POLDER, MISR, or Along Track Scanning Radiometer 2 
(ATSR-2) characterize aerosols by using the different look angles to measure surface 
reflectance and use lookup tables to model the atmosphere (Veefkind et al. 1998, 
Martonchik et al. 2004, Deuze et al. 2000). For example, the dual view algorithm, 
developed by Veefkind et al. (1998), uses two look angles to separate the atmospheric 
and surface contributions to the satellite radiance. The multi-angle method, though quite 
accurate, includes a step in the algorithm to identify and ignore regions with high surface 
reflectance (Robles Gonzalez et al. 1997). 

B. THE SHADOW METHOD 

Except where otherwise cited, the information provided in this description of the 
shadow method was taken from Vincent (2006). The shadow method exploits high- 
resolution imagery to measure the radiance difference between building shadows and 
adjacent sunlit regions to calculate total optical depth. Several commercially-available 
satellites have sufficient horizontal resolution, from 46 cm for the WorldView-1 
panchromatic channel to 2.44 m for the QuickBird multispectral channels (DigitalGlobe, 
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cited June 2010a/b), to use the shadow method. Figure 2 depicts the basic components of 
radiance received at a satellite, from a shadow and an adjacent sunlit region where both 
occur over an identical surface. Direct transmission, diffuse transmission, and diffuse 
reflection components make up the sunlit signal, while the shadow signal lacks the direct 
transmission component. For a homogeneous surface, where a portion is covered by 
shadow, subtracting the shadow radiance from the sunlit radiance isolates the direct 
transmission component of the signal, the portion of the signal that is the key to 
calculating optical depth. 



Figure 2. Depiction of the difference between the shadow and sunlit signal received by 
a high-resolution satellite. Subtracting the two signals isolates the direct 
transmission component of the signal, which can be used to calculate total optical 

depth (from Vincent 2006). 

1. General Calculation 

The equation for calculating TOD is derived from Liou’s (2002) development of 
the principles-of-invariance (POI) method to describe the multi-scattering problem. 
Vincent (2006) provides a detailed derivation of the governing equation: 
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where d> is the single scatter albedo, P is the scattering phase function (represented by 
the Henyey-Greenstein phase function), g is the asymmetry parameter (0.65), ©is the 
scattering angle, and (p',^')and (p",^")are the cosine of the incident and resultant zenith 


and azimuth angles (Henyey and Greenstein 1941). The single scatter albedo value 
varies regionally, but Vincent used values of 0.88 for Beijing, China and 0.95 for other 
locations. Additionally, Vincent set the asymmetry parameter value to 0.65 to account 
for both aerosol (0.6-0.7) and Rayleigh particles (~ 0). 


Y approximates the diffuse sky radiance. Although the diffuse field consists of 
aerosol and molecular affects, the molecular affects were ignored when establishing the 
asymmetry parameter in Equation (3) because aerosol affects dominate at all but the 
lowest optical depths. Had the asymmetry parameter included the molecular affects, its 
value would decrease. Vincent concluded this would cause the scattering phase function 
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and MAR to increase, thereby lowering the AOD and increasing the low bias of the 
shadow method. Therefore, molecular affects were only accounted for during the final 
partitioning of the total optical depth. 

To estimate AOD, the contribution from Rayleigh scattering, or molecular 
Rayleigh optical depth (6 R ), must be calculated and subtracted from the TOD. The effect 
of molecular absorption is considered negligible since most satellite sensor wavelength 
bands operate in atmospheric transmission windows. The molecules responsible for 
Rayleigh scattering are standard atmospheric constituents and may be calculated globally 
for a given wavelength using equations defined by Russell et al. (1993): 

5 r (A,) = (0.00864 + 6.5 x10~ 6 h) U bW (p / p 0 ) (5) 

b(k) = 3.916 + 0.074k + 0.050 /k (6) 

where H is the height above sea level of the radiometer in kilometers, p is the 
atmospheric pressure at H in hPa, and po is the sea-level reference pressure (1013.25 
hPa). The depolarization factor (0.00864) was chosen based on the best fit data (Frolich 
and Shaw 1980, Young 1980). Using Equations (5) and (6), 5r was approximated using 
the center effective wavelength of the satellite channel. 

2. Sensitivities 

Vincent’s sensitivity analysis of key terms in Equation (1) demonstrated the 
shadow method is most sensitive to the radiance difference between sunlit and shadow 
regions and surface reflectance and much less sensitive to 7, solar zenith angle, and 
sensor zenith angle. Additionally, Vincent detailed the expected impact on the shadow 
method when cumuliform, stratiform, or cirrus clouds were present. In the case of 
cumuliform or stratifonn clouds, images with cloud coverage of more than 20% were 
avoided. Thin cirrus was more difficult to detect using visible channels and resulted in 
retrieved radiances at the sensor being higher, causing an overall high bias to the 
calculated AOD. 
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3. 


Limitations and Target Selection Criteria 


Vincent also found several limitations to the shadow technique. First, the 
minimum allowed radiance difference was 10 W m" sr' nm' (RU) and the minimum 
allowed surface reflectance was 0.15. These thresholds were set because of the high 
sensitivity to changes in radiance difference and surface reflectance below those values. 
Such sensitivity leads to large errors in AOD calculations with small errors in radiance 
difference or surface reflectance. The range of calculable AOD was set from 0.1 to 2.0. 
The minimum value was established to avoid negative AOD after subtracting 5r 
(Rivenbark 2009). The maximum was detennined by calculating TOD using the 
combination of allowable values that gave the maximum result. Furthennore, when 
optical depth is above 2.0 the surface is obscured and shadows are indistinguishable. 

The method also required careful target selection in urban environments and only 
blue and NIR channels provided reliable results over grass backgrounds. This resulted in 
several recommended target selection criteria. These criteria included using large 
shadows over a homogeneous surface with a minimum shadow size of 4 pixels wide and 
4 pixels deep (from about 1 to 32 square meters, depending on satellite resolution), using 
shadows that are at least 1 km away from other structures, and avoiding highly reflective 
shadowing structures and the subsequent adjacency effects. Vincent suggested two 
approaches to automating the technique, dynamic and fixed. The dynamic approach 
relied on isolating shadows of opportunity while the fixed approach required the 
preidentification of buildings meeting specified criteria for analysis. 

4. Other Observations 

Regardless of channel or background, the shadow technique was a manually 

intensive procedure with more time being spent identifying good shadows than 

performing the actual calculation. Process streamlining perfonned by Rivenbark (2009) 

and Dombrock (2007) took steps to improve the calculation time while exploring aspects 

of shadows that contributed to error in optical depth. In each study, errors in the code 

written in Interactive Data Language (IDL) impacted their findings. For example, when 

calculating the scattering phase function as part of the Y calculation in Equation (3), a 
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syntax error in code raised the denominator to the power of 1 instead of 1.5. The overall 
impact of this mistake is difficult to determine mathematically, since the error can cause 
either an increase or a decrease in the value of the scattering phase function, depending 
on the scattering angle. Also, additional iterations in the azimuth angle loop nests in 
Equation (2) contributed to the low bias by causing higher values of T and lower values 
of optical depth. 

C. SHADOW DETECTION IN SATELLITE IMAGERY 


The key step in developing a fully automated approach to the shadow method is 
finding suitable shadows within the image. One satellite-based method of finding 
shadows uses color space indices (Sarabandi et al. 2004, Gevers and Smeulders 1999). 
The color space indices are defined as follows: 

_ f R 

C, = arctan - 

v max(G,B) 

_ f G ) 

C, = arctan - 

v max(R,B) y 


C, = arctan - 

l v max(R,G) 

where R, G, and B represent the radiance from the red, green and blue color bands, 
respectively. Coupling the output from these indices with a 3x3 texture filter clearly 
outlines the boundary between a shadow and sunlit region. Sarabandi et al. (2004) 
determined that C 3 worked best for locating shadows. Due to more Rayleigh scattering 
of blue light than red or green light, C3 is larger in shadow regions than sunlit regions. 

Although capable of identifying shadows, this satellite-based approach is poorly 
suited for automating the shadow method for two reasons. First, C3 may be less effective 
in regions with large aerosol concentrations since the amount of scattered red and green 
light will increase. Second, this technique lacks necessary infonnation to select sunlit 
pixels with the same surface reflectance as shadow pixels, particularly when using 
panchromatic channels, where surface reflectance estimation techniques are not available. 
Similar problems would likely arise with other satellite-based approaches. The next 
section describes the FASTEC procedures for identifying usable shadows and calculating 
AOD. 
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III. METHODOLOGY 


An automated approach, based on trigonometry, is created using the shadow 
technique and overcoming the shortfalls of satellite-based shadow identification 
techniques mentioned in Chapter II, Section C. Knowing pixel elevations and the 
azimuth and elevation of the sun and satellite, the process for identifying shadow and 
sunlit pixels is relatively straightforward. If the line-of-sight between the sun and a pixel 
(called the target pixel) is blocked by a higher elevation pixel, a target pixel is in shadow. 
If not, a target pixel is sunlit (see Figure 4). Similarly, the new approach ensures the 
satellite does not collect radiance off an obstruction, instead of the target by ensuring the 
line-of-sight between the satellite and the target pixel is not blocked. 

A. TARGET SELECTION CRITERIA TESTING 

Past research identified several target selection criteria for identifying good 
shadow candidates for the shadow method (Vincent 2006, Dombrock 2007, Rivenbark 
2009). Criteria include avoiding partial shadow pixels along the shadow edge and 
avoiding pixels where surrounding objects limit diffuse sky radiance. Detailed 
descriptions of how the automated approach addresses these concerns are described in 
Chapter III, Section C. The hypotheses listed below describe key target selection criteria. 

Before these criteria were tested, variables used to discriminate between shadow 
generators were tested to eliminate tree shadows from the sample. Tree shadows may 
allow some sunlight to penetrate through the leaf canopy into the shadow, increasing the 
shadow radiance, invalidating the theory behind the shadow technique (see Figure 2). The 
first of these variables sets the minimum height of the shadow generator and the second 
sets the minimum number of pixels required for AOD calculation. To determine the 
combination of these variables that limits the sample to buildings, different values of the 
minimum shadow generator height were applied to the GSFC tile and the output 
examined to find the value that leaves only building shadows. Furthermore, plots of 
AOD for the entire DSM were created to show how it varies as tree shadows are 
introduced to the sample. 
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Figure 3. Example of pixel rejections from FASTEC. Orange represents the building, 

black is the shadow, and gray is rejected shadow. Panel 1 is the full shadow, 2 is 
after edge rejection, and 3 is after building proximity rejection. 

1. Hypothesis I: It Is Important to Avoid Partial Shadow Pixels Along 
the Shadow Edge 

Subjective examination of satellite imagery shows that radiance values transition 
across several pixels along the shadow edge. In the transition region, shadow radiance is 
slightly higher and sunlit radiance is slightly lower than pixels that are fully in shadow or 
sunlit. Including these pixels in the optical depth calculation, lowers the radiance 
difference and increases AOD. The edge pixel rejection variable enables the automated 
technique to eliminate partial shadow pixels along the shadow edge (see Figure 3). To 
test this hypothesis, plots of radiance based on distance from shadow edge were reviewed 
to determine the extent of the influence of the shadow edge on the shadow region. 
Additionally, the edge pixel rejection threshold was varied and the resulting AOD values 
compared for individual shadows and an entire DSM file. Chapter III, Section B.5.C. 
describes the algorithm FASTEC employs to eliminate partial shadow edge pixels. 

2. Hypothesis II: It Is Important to Choose Shadow Pixels Far from the 
Shadow Generator to Calculate Optical Depth 

Based on the past research, it is suggested the shadow region should only include 
pixels far from the shadow generator. In theory, reflection off surrounding objects may 
increase pixel radiance while blocking by surrounding objects may decrease pixel 
radiance. Regardless of which modification dominates, the total signal is no longer 
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composed of the basic components listed in Figure 2. By rejecting a percentage of the 
shadow pixels near the shadow generator based on the size of the overall shadow (see 
Figure 3) this adjacency affect is reduced. To test this hypothesis, the correlation 
between the pixel radiance and distance from the shadow generator was calculated for 
several shadows. Also, the building percentage rejection threshold was varied on 
individual shadows and on an entire DSM file to examine the impact on AOD. Chapter 
III, Section B.5.a. describes the algorithm FASTEC employs to eliminate pixels too close 
to the shadow generator. 

3. Hypothesis III: Sky Dome Blockage Is Important When Choosing 
Valid Shadow and Sunlit Pixels for Calculating Optical Depth 

Based on the same reasons as described for Hypothesis II, selecting pixels based 
on sky dome blockage may be important to the shadow technique. In order to determine 
the influence of sky dome blockage on the radiance of a pixel, the correlation between 
them was calculated for several shadows. Additionally, AOD was calculated for different 
maximum sky dome blockage thresholds for individual shadows and an entire DSM file. 
Chapter III, Section B.5.b. describes the algorithm FASTEC employs to eliminate pixels 
based on sky dome blockage. 

4. Hypothesis IV: Eliminating the Extreme Radiance Values from the 
Sample Improves AOD Calculations 

Radiance variation exists in the shadow and sunlit regions due to variations in 
surface reflectance and position differences between the satellite and DSM files. Surface 
reflectance variations may be caused by different surfaces or different objects such as 
trees, sidewalks, or cars. Position differences between the satellite and DSM files are 
partially corrected by the offset correction algorithm, but some error in position, rotation, 
or distortion may still exist. Therefore, FASTEC eliminates percentages of the extreme 
radiance values within shadow and sunlit regions to try and ensure the sampled regions 
are representative of the true shadow/sunlit radiances. Assuming FASTEC correctly 
identifies the shadow, the bulk of the pixels within the samples are representative of the 
true shadow/sunlit radiances. Sorted values of pixel radiance were plotted to demonstrate 
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this problem. Also, thresholds for eliminating extreme values of radiance from both the 
shadow and sunlit region were varied with resulting AOD values plotted for both 
individual shadows and across an entire DSM. Chapter III, Section B.5.d. describes the 
method used to eliminate extreme radiance pixels. 

5. Hypothesis V: It Is Possible to Get an Idea of Uncertainty When 
Calculating AOD on Entire Scenes. 

After setting target selection criteria based on the tests described by Hypotheses I- 
IV, calculations on an entire scene will provide a number of shadow/sunlit pairs and 
computed AOD values. The range of AOD values in the scene will give an idea of the 
uncertainty of the automated technique. To test, AOD values were calculated for all 
available imagery using entire DSM files. In cases where ground truth data is also 
available, AOD values for each shadow were plotted against ground truth and bias was 
calculated. The data is analyzed by region (for each DSM file) as well as across all 
regions (all DSM files). In the case of multispectral imagery, plots of AOD between 
channels were examined. Since the overlap between most of the available imagery and 
the DSM files contains several qualified shadow generators, a range of AOD values was 
available for most test cases. 

6. Hypothesis VI: The DSM Can Have Significantly Lower Resolution 
Than the Satellite Data and Still Work Effectively 

The computational expense of the automated algorithm is most influenced by the 
size of the DSM file, limiting the size of the region over which the automated method can 
be run operationally. Using lower resolution imagery allows calculations over a larger 
area without increasing computation time. For example, increasing the horizontal 
resolution of the DSM file from one meter to two meters, quadruples the areal coverage 
of the DSM with the same computational cost. As long as the resolution is sufficient to 
identify shadow regions with several pixels, it is believed the resulting AOD values will 
be comparable to the full resolution DSM. To test this hypothesis, results with varying 
horizontal resolution DSM files were compared to the full resolution results. 
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B. TRIGONOMETRIC SHADOW DETECTION METHOD 



Figure 4. Trigonometric approach employed by FASTEC to find shadow pixels based 
on elevation. Red arrows indicate shadow pixels whose path to the sun is blocked 
by a higher elevation pixel. The green arrow indicates an unblocked sunlit pixel. 
Black pixels are obstructed from view of the satellite by the building. 


A high-resolution surface elevation map, known as a Digital Surface Model 
(DSM), provides the critical data for FASTEC. Using this data instead of satellite data, 
ensures that the automated approach works in the presence of obscurations such as clouds 
and aerosols. To execute, FASTEC requires overlapping DSM and satellite files as 
inputs. Output includes a list of shadows with total optical depth (TOD) and AOD values 
for each qualifying shadow in the DSM file and general details of shadow/sunlit selection 
criteria. Figure 5 depicts a top-level flowchart of the automation program. Subroutines 
for each block in this flowchart are described in detail in the rest of this section. 
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Figure 5. Top-level flowchart of FASTEC. 

1. Preparation 

The preparation stage is made up of two steps with step one being performed by 
the user prior to running the program and step two being performed by the automation 
program during execution. Before the automated method can begin, the DSM and 
satellite files are converted into flat file format, composed of a header file and a data file. 
This makes it easy to retrieve key header information. The details of this step as well as a 
description of the input files are detailed in Chapter III, Section C. 1 and 2. 

The automated portion of the preparation stage begins by reading in a host of 
variables from a configuration file. The configuration file allows the user to quickly 
change variable values and eliminates the need to recompile the program whenever these 
variable values are changed. For this program, the configuration file variables include the 
input and output file names, shadow and sunlit modification variables, offset correction 
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and dome blockage criteria, and AOD calculation variables (see Appendix A for a 
description of each variable). Next, key information from the DSM and satellite flat file 
headers are read into the program including: the number of rows and columns of data, 
the number of channels of data, the data storage format, the latitude and longitude of the 
top left pixel, and pixel size in both the north-south and west-east directions. 
Additionally, necessary infonnation from the image metadata file is read including: the 
absolute calibration factor and effective bandwidth of each channel, the date and time of 
the image, the satellite name, and the sun and satellite azimuths and elevations. Equation 
(8) converts Digital Numbers (DN) into spectral radiance using absolute calibration 
factor and effective bandwidth: 


L,(i) = DN, 


absCalFactj 

effBandwidthj 


( 8 ) 


where Eli) is the spectral radiance, i is the channel number, and absCalFactj and 


effBandwidthj are the absolute calibration factor and effective bandwidth, respectively. 
DN is essentially the relative spectral radiance and is based on the amount of energy 
detected by the sensor. Additionally, the values for §r and spectral solar irradiance are 
chosen from Table 1 during this step. § R values were calculated by taking a response- 
weighted average using Equations (5), (6), and (9): 


8 rt — 


j*8 R (A.)0(A.)dA. ^(k) 


(9) 


where 5 RT is the total MROD across the channel and cp(A) is the specific spectral 


response for that wavelength and satellite channel. An example of the spectral response 
curve for the QuickBird satellite is provided in Figure 6. The summation in Equation (9) 
was performed using a 50 mn wavelength interval across the channel. Spectral solar 
irradiance was calculated using the same approach where the wavelength-specific 
spectral solar irradiances were taken from Wehrli (1985). 
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Sat 

Band 

Wavelength 
Range (nm) 

Center Effective 
Wavelength (nm) 

Molecular Rayleigh 
Optical Depth 

Spectral Solar 
Irradiance 
fW m ' 2 nm' 1 ) 

Quickbird 

Blue 

450-520 

486 

0.172 

1940 

Green 

520-600 

547 

0.105 

1851 

Red 

630-690 

652 

0.051 

1593 

Near-Infrared 

760-900 

804 

0.026 

1127 

Panchromatic 

445-900 

689 

0.061 

1558 

Ikonos 

Blue 

445-516 

497 

0.165 

1916 

Green 

506-595 

560 

0.097 

1835 

Red 

632-698 

666 

0.047 

1539 

Near-Infrared 

757-853 

792 

0.027 

1170 

Panchromatic 

526-929 

727 

0.048 

1468 

WV1 

Panchromatic 

400-900 

667 

0.070 

1587 


Table 1. List of Molecular Rayleigh Optical Depth and Spectral irradiance 

constants by satellite and channel. 


The last steps in the preparation stage expand the row and column sub-region (if a 
sub-region was defined) and reads in the entire DSM file. Expanding the sub-region 
enables FASTEC to identify usable sunlit pixels outside the user-defined area to pair with 
shadows inside the sub-region. The entire DSM file data is stored in an array for the step 
described in the next section. 




Wavelength 


Figure 6. Relative spectral response for the QuickBird satellite derived from data 

provided by DigitalGlobe. 
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2. 


Identify Shadow and Sunlit Pixels 


This step identifies shadow and sunlit pixels and determines which pixels are 
blocked from the view of the satellite. The results are stored in an array and used by 
subsequent steps in the program. The results array contains one of three possible values: 
distance from the shadow generator for in-shadow pixels, zero for sunlit pixels, or invalid 
for pixels blocked from view by the satellite. 



Figure 7. Calculation of in-shadow pixel h(r t ar,Ct ar ). The height of the beam (hb) 
between the target pixel and the sun is less than the height of pixel h(ri os ,ci os ). Diff 
is the difference between h(ri os ,ci os ) and hb. 


Pixels are considered shadow whenever the line-of-sight between the target pixel 
and the sun is blocked by a higher elevation pixel (see Figure 7). The elevation of the 
line-of-sight between the sun and target pixel is computed using the below equations: 


dist = J[dy(r l0S -r tar )]'+[dx( 


C los C tar 


h b = h (Far > Car ) + dist x tan (elev sun ) 


( 10 ) 

( 11 ) 


where hb is the elevation at pixel (ri os , ci os ) of the line-of-sight between the sun and the 
target pixel (r ta r,Ctar), dx (dy) is the pixel size in meters in the east/west (north/south) 
direction, h(r tar ,c tar ) is the height of the target pixel, and elev sun is the elevation angle of 
the sun in radians. This calculation is repeated for every pixel on the azimuth between 
the target pixel and the sun until hb exceeds the elevation of the tallest pixel in the DSM. 
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Shadow pixels store the distance from the target to the shadow generator in the 
result array. Some shadow pixels, such as those having two tall buildings between them 
and the sun, may have multiple shadow generators. In that situation, the result array 
value is determined by the pixel with the largest difference (diff) between the shadow 
generator height h(ri os ,ci os ) and h b (see Figure 7). This ensures the stored value coincides 
with the shadow generator that places the target pixel deepest in shadow. This approach 
will cause an occasional good shadow pixel to be marked invalid when eliminating 
unwanted pixels (see Chapter III, Section B.5.a). Any pixel whose line-of-sight is 
unblocked is considered sunlit and is labeled with a zero in the result array. 

After identifying each pixel in the DSM as shadow or sunlit, the view between the 
satellite and the ground is considered. This ensures that pixels representing objects 
obstructing the surface when viewed from the satellite are identified and labeled invalid 
in the result array. The procedure for perfonning this test uses a similar approach as 
described above with the satellite location substituted for that of the sun. Aside from 
that, the only difference is that any obstructed surface region is immediately marked 
invalid in the result array. When this step completes, the result array contains pixels 
identified by distance from shadow generator for shadow pixels, zero for sunlit pixels, or 
invalid for pixels obstructed from view of the satellite. 

3. Identify Coherent Shadows 

Once shadow and sunlit pixels are identified, FASTEC groups adjacent shadow 
pixels as one shadow and assigns that shadow a number, sequentially from one, for 
identification. The procedure begins by examining the result array from the previous 
section (any pixel whose value is greater than zero). Once a shadow pixel is found a 
recursive subroutine is called to find neighboring shadow pixels. Any neighboring 
shadow pixel that is found is identified as belonging to the same shadow and the same 
recursive subroutine is repeated for that shadow pixel. This process is repeated until 
every shadow pixel has been assigned a shadow number. To improve perfonnance in 
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later steps, the range of rows and columns and the pixel count for each coherent shadow 
are also stored. The coherent shadow array contains positive integers for shadow pixels, 
zero for sunlit pixels, and invalid for other pixels. 

4. Correct for Pixel Offset Between DSM and Satellite File 

The latitude and longitude of each point in the DSM and satellite file must match 
precisely for FASTEC to properly measure shadow and sunlit radiance and calculate 
AOD. Due to the poor geolocation of the satellite imagery used in this study, pixel 
locations within these files rarely match and correction is required. To demonstrate the 
impact of this problem, shadow pixels identified by the automated method were 
compared to those identified by hand in a satellite file (see Figure 8). In Figure 8, panel 
A represents the results from FASTEC, where green identifies shadow pixels and yellow 
identifies partial shadow pixels. Due to resolution differences between the DSM file (one 
meter) and the satellite file (2.4 meters), the resolution from the automated results in 
Panel A was matched to the satellite resolution in Panel B. This led to between four and 
nine DSM pixels comprising each pixel in Panel A. In cases where every DSM pixel 
within a satellite pixel was identified as shadow by FASTEC, the comparison identified 
the pixel as fully in shadow. In cases with a mixture of shadow and sunlit pixels existed, 
the comparison identified them as partial shadow. Panel B represents the results of a by¬ 
hand analysis of the satellite image. Panel C shows the shadow pixel overlap (162 pixels) 
between Panels A and B without performing offset correction. By shifting the latitude 
and longitude of every pixel in the DSM file equally, the overlap increased to 218 pixels, 
as shown in Panel D. Without offset correction, 38% of the FASTEC shadow pixels 
were actually sunlit, resulting in higher shadow radiance and higher values of optical 
depth. This case required an eight meter shift, performed by hand, to maximize the 
overlap. Other cases required shifts of 100 meters or more, most likely due to the 
imprecise orthorectification applied to the satellite files. 
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Figure 8. Four-panel plot of offset between satellite and DSM files for Goddard Space 
Flight Center on October 26, 2006. In Panels A and B, green is shadow and 
yellow is partial shadow pixels. Panel A displays the results from FASTEC and B 
displays the results from a by-hand analysis of the satellite file. Panels C displays 
overlap between A and B without offset correction (162 pixels). Panel D displays 
overlap with offset correction (218 pixels). 

Regardless of the amount of shift required, correcting for the offset between the 
input files is critical to FASTEC. An automated offset correction procedure that adjusts 
the latitude and longitude of the upper left corner in the DSM file to find the lowest 
average shadow pixel radiance was developed. The range over which the average 
radiances are computed is set in the configuration file. Shadow pixels are identified 
using the coherent shadow array, with correction for edge contamination (see Chapter III, 
Section B.5.b.). This procedure corrects for translation error between the DSM and 
satellite file. It does not correct for rotation error or distortion error. To improve 
computation time, the correction is performed with a coarse and fine stage. During 
coarse correction, the coordinates of the upper left corner of the DSM are shifted by an 
increment in the configuration file, nominally 4 pixels in both the north/south and 
east/west directions. After each shift, the average coherent shadow radiance of one 
satellite channel is calculated for each valid target shadow. The offset that provided the 
lowest average radiance is the starting point for fine correction. During fine correction, 
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the process is repeated for every possible shift between the coarse increment and the 
coarse correction best fit. After the lowest average shadow radiance is found, the latitude 
and longitude of the DSM file is adjusted by that offset interval. 

5. Eliminate Unwanted Pixels 

Several procedures were developed to eliminate pixels whose radiances are not 
representative of the radiance in shadow or sunlit regions because of reduced diffuse 
radiance, combined shadow/sunlit radiance, and differing surface reflectance. Each 
procedure is geared toward testing target selection theories as outlined in Vincent (2006), 
Rivenbark (2009), and Dombrock (2007). Figure 3 displayed examples of output from 
FASTEC for two of these criteria. Each of the rejection methods are described in more 
detail in the remainder of this section. Different variables in the configuration file control 
execution of these procedures, thereby making them easily adjustable. 

a. Pixels That Are Too Close to the Building 

As portions of the sky are blocked by taller structures, the amount of 
diffuse sky radiance reflecting off a pixel is reduced, thereby lowering radiance values 
and increasing error in AOD measurements. Rivenbark (2009) and Dombrock (2007) 
suggested selecting shadow pixels far from the shadow generator in order to mitigate the 
impact of this effect. FASTEC includes two methods of eliminating pixels based on their 
proximity to the shadow generator; the first eliminates pixels using their physical distance 
from the shadow generator while the second eliminates a percentage of the pixels within 
the overall shadow (see Figure 9). 
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distp* 


Figure 9. Identifying pixels by proximity to shadow generator. h b i dg is the height of the 
building, h b is the height of the line-of-sight between the target pixel and sun, 
distpi x is the distance from the shadow generator to the target pixel, dist max is the 
farthest extent of shadow, and p user is the percentage of pixels to mark invalid. 

In the interest of computational efficiency, these calculations are actually 
performed when identifying shadow and sunlit pixels (see Chapter III, Section B.2) but 
are stored in a separate array. These adjustments are separated to prevent error in the 
offset correction algorithm. Additionally, it allows the output from this step to be stored 
in a file and examined separately during testing. Assuming the radiance is lowest near 
the building where diffuse sky radiance is blocked by the building, performing the offset 
correction after this step would choose offset values that place the shadow region too 
close to the building, the exact location this step is trying to avoid. This would cause the 
alignment between the DSM and satellite file to be incorrect, likely resulting in 
significantly lower sunlit radiance values as many pixels identified as sunlit in the DSM 
would fall over shadow in the satellite file. 

When selecting pixels based on the physical distance from the shadow 
generator, the distance from the target pixel to the shadow generator is compared to a 
threshold as defined in the configuration file. If the distance is less than the threshold, the 
shadow pixel is marked invalid. To select pixels based on their proximity by percentage 
of the overall shadow, the ratio displayed in Figure 9 is compared against the building 
percentage rejection threshold in the configuration file. If the ratio is less than the 
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threshold, the pixel is marked invalid. The results from this step are integrated into the 
coherent shadow array immediately after offset correction is performed. 


b. Pixels Whose Sky Dome Is Blocked by Adjacent Pixels 


As described in the preceding section, diffuse sky radiance is blocked by 
any tall objects near a pixel. The test based on distance from the building addressed the 
problem between shadow pixels and the shadow generator. To fully address the impact 
of diffuse sky radiance, the automated approach must examine other surrounding objects. 
A procedure was developed to examine the sky dome blockage of shadow and sunlit 
pixels before allowing them to be used in the AOD calculation. The sky dome is 
considered the region of the sky that is visible to the target pixel. 



Figure 10. Celestial dome blockage. The dark line depicts the angle representing dome 

blockage along that azimuth. 


The procedure determines the percentage of sky dome blockage by 
averaging the blockage along equally-spaced azimuths (see Figure 10). The angle 
between every pixel within the search range and the target pixel is computed along each 
azimuth in a circle around the target pixel. The azimuth interval, distance interval, and 
search range are assigned using the configuration file. The angle (9 ) between the pixels 
is computed using the following equation: 

f \ 


9 = tan 1 


h(r 2 ,c 2 )-h(r tar ,c t J 
■x/dyif -fa,-) 2 +dx(c 2 -c lar ) 2 


180 


( 12 ) 
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where (r 2 ,C2) and (r ta r,Ctar) are the along-azimuth pixel and target pixel respectively, h(r,c) 
is the height of that point, and dx and dy are the pixel resolution in the north/south and 
east/west directions. The average maximum blockage angle along each azimuth is 
divided by 90 to get the overall blocked percentage of the sky dome. If the average 
blockage exceeds the threshold in the configuration file, the pixel in the coherent shadow 
array is labeled invalid. The output from this procedure is integrated into the coherent 
shadow array prior to offset correction. Additionally, the average blockage calculation 
results are stored in an “average blockage” array. 

c. Pixels That Are Too Close to the Shadow Edge 

Pixels along shadow edges are partially sunlit. This is apparent because 
the radiance values increase significantly near the edge of the shadow. Figure 11 
displays the shadow pixel radiances as detennined by FASTEC. In this example, pixels 
near the center of the shadow have spectral radiance values between 15 and 17 RU while 
radiances get larger as the shadow edge is neared. This pattern is reversed for sunlit 
pixels. Including these edge pixels raises shadow radiance and lowers sunlit radiance, 
causing computed AOD values to increase. In the satellite image, one would expect this 
contamination to affect only those pixels along the shadow edge. Due to the nearest 
neighbor method of retrieving radiances for each pixel and issues with offset between the 
images, this problem may impact more than just those pixels closest to the shadow edge. 
A procedure was created to eliminate contaminated edge pixels from the computation. 
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Figure 11. Example of shadow radiances from FASTEC. Spectral radiance values are 

e 2 11 

units of W m" sr‘ nm~ (RU). Building location is estimated. 


The procedure to eliminate shadow edge pixels examines every pixel 
within range (as assigned by the configuration file) of valid shadow pixels. If a sunlit 


(shadow) pixel is within range of the target shadow (sunlit) pixel, the target pixel is 
labeled invalid. The result from this step is stored in an “interior” array for later use. 
The interior array is integrated into the coherent shadow array prior to offset correction to 
reduce the shadows used in offset correction to those most likely to have optical depth 


calculations performed. 


d. Pixels with Elevation Issues or Extreme Radiances 


The shadow method requires the shadow and sunlit regions to have the 
same surface reflectance. Without surface reflectance maps of comparable resolution. 


this requirement is impossible to fully satisfy. In order to improve the likelihood this 
requirement is met, two final corrections are performed. 
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The first correction examines the elevation of every point within the 
shadow region. If the difference between the maximum and minimum elevation exceeds 
a threshold in the configuration file, the pixel whose elevation falls farthest from the 
mean elevation of the shadow region is labeled invalid. This process is repeated until the 
difference between the max and min elevation no longer exceeds the threshold. This 
procedure eliminates trees, small buildings, cars, or other tall objects that most likely 
have different surface reflectance within the shadow region. It does not account for 
different surface types, such as asphalt and grass, which change surface reflectance. A 
similar procedure for sunlit pixels is described in the next section. 

The second correction is the only correction that is based on satellite pixel 
radiance. It assigns radiances to the remaining valid shadow pixels, sorts them by 
radiance, and then eliminates percentages of the extremes based on thresholds in the 
configuration file (see Figure 12). The pixels are sorted by radiance as they are read into 
the array and the extremes are not used to compute the average radiance of the shadow. 
This procedure eliminates small areas in the shadow region with different surface 
reflectance and counteracts some offset correction error. 



Figure 12. Histogram of radiance values with extreme radiance rejection. The actual 
distribution is not necessarily normal as depicted in this graphic. 
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6. 


Identify Accompanying Sunlit Pixels 


For each shadow that is large enough (as defined in the configuration file) for 
AOD computation, sunlit radiance is required. The procedure to compute sunlit radiance 
begins by associating a set of candidate sunlit pixels based on their proximity to the 
shadow. The search region for valid sunlit pixels is limited to a 180 degree region around 
the shadow, perpendicular to the sun (see Figure 13). Using an azimuth and distance 
approach similar to the procedure described in the section for sky dome blockage, this 
procedure searches for candidate sunlit pixels. 


Non-s hadow R egion 

Shadow 

pixels 


Figure 13. Region for sunlit pixel search. 

After identifying the initial sunlit candidates, the procedure eliminates unwanted 
pixels. First, sunlit pixels are eliminated based on their sky dome blockage, using the 
same procedure described in Chapter III, Section B.5.b. Next, it eliminates pixels based 
on elevation using a similar process to the method described in the previous section. 
Unlike the shadow pixel elevation elimination, this process compares the elevation of 
each sunlit candidate to the average elevation of the shadow region. Any pixel elevation 
that falls outside the tolerance, as defined in the configuration file, is eliminated from the 
sunlit pixels dataset. Since it is possible to associate a sunlit pixel with more than one 
shadow, this elimination is done for each individual shadow. This approach is done to 
align the shadow and sunlit regions in an effort to improve confidence that they have the 
same surface reflectance. 
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The final elimination step removes the extremes of the sorted radiances using the 
same approach described in Chapter III, Section B.5.d. Finally, the average sunlit 
radiance is computed from the remaining sunlit pixels. The automated shadow detection 
program described above greatly increases the total number of shadow/sunlit pairs for a 
given time period when compared to manual calculations. The last major section of the 
program ties in past research to this automated shadow detection algorithm. 

7. Calculate AOD for Each Shadow 

The AOD calculation is performed for each shadow/sunlit pair that meet the 
minimum size requirement using a series of subroutines to compute Equations (1) - (4). 
The subroutines were optimized to improve computation speed but follow the 
methodology used in past research (Vincent 2006, Dombrock 2007, Rivenbark 2009) 
with corrections to the issues described in Chapter II, Section B.4. In instances where the 
shadow radiance is larger than the sunlit radiance (due to poor target selection, offset 
correction, or cloud cover), the calculation is not performed and dummy TOD and AOD 
values are passed to the output file to allow for further examination of those shadows. 
The output from the automated algorithm is described in detail in Chapter III, Section 
C.3. 

C. DATA 


DSM-LIDAR Satellite 



Figure 14. Examples of input data used in this research. Both images are over the 
Goddard Space Flight Center (GSFC) near Washington D.C. The QuickBird 
satellite image is from October 26, 2006. 
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1 . 


Digital Surface Model (DSM) 


For the purposes of this study, the term Digital Surface Model (DSM) refers to the 
one meter resolution elevation maps used as the basis of the automation technique. In 
these maps, each data point (or pixel) contains the elevation of the first solid object 
encountered when looking down at the surface from above. Several methods have been 
developed for creating a DSM. Two measures are commonly employed to describe 
accuracy; Circular Error of 90% (CE90) for horizontal accuracy and Linear Error of 90% 
(LE90) for vertical accuracy. CE90 is the minimum diameter of the horizontal circle that 
can be centered on all photo-identifiable Ground Control Points (GCPs) and also contain 
90% of their respective twin counterparts acquired in an independent geodetic survey. 
LE90 is the linear vertical distance that 90% of control points and their respective twin 
matching counterparts acquired in an independent geodetic survey should be found from 
each other (GeoVAR 2009). Two techniques for creating high resolution DSM data are 
stereo-pairs and light detection and ranging (LIDAR). 

Stereo pairs consist of two images of the same scene taken at near simultaneous 
times with slightly different view angles. Ground control points provide baseline 
location and elevation infonnation to calibrate the image and improve its accuracy. 
DigitalGlobe advertises accuracy of 6.5m CE90 and 6.5m LE90 (DigitalGlobe, cited June 
2010c) although Low et al. (2006) improved CE90 to 5.0m or better. 

Another method of creating a DSM uses LIDAR instruments flown aboard 
aircraft. The U.S. Army Topographic Engineering Center (TEC) BuckEye system 
employs a LIDAR sensor to create one meter resolution DSM data of 0.5m CE90 and 
0.3m LE90 accuracy (TEC 2009). The high resolution and accuracy of this data makes it 
a superior choice for the automation technique. For this study, the LIDAR data is 
distributed as approximately 5km x 5km tiles in GeoTiff file format. ENVI 4.5 is 
visualization software, capable of displaying and manipulating data in a multitude of 
formats. This software is used to convert these DSM files into flat files to simplify data 
ingest. 
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2 . 


Satellite Data 


High-resolution commercial satellite data from the IKONOS, QuickBird and 
WorldView-1 satellites were used in this study. QuickBird offers panchromatic and 
multispectral imagery with CE90 of 23m and horizontal resolutions of 60cm and 2.4m, 
respectively. WorldView-1 offers a panchromatic imagery with CE90 of 6.5m and 
horizontal resolution of approximately 50cm (DigitalGlobe, cited June 2010a/b). 
IKONOS offers multispectral and panchromatic imagery with CE90 of about 6m and 
horizontal resolutions of lm and 3.28m, respectively (Dial et al. 2003, IKONOS, cited 
May 2010). All satellite data was ordered via the Internet using the Web-based access 
and retrieval portal on the National Geospatial-Intelligence Agency Web site in NITF 2.0 
format. This study used basic imagery from DigitalGlobe (QuickBird, WorldView-1). 
This imagery is not orthorecti lied although a built-in process in ENVI, using rational 
polynomial coefficients, was used in an effort to remove some distortion. ENVI software 
was also used to save the imagery as flat files. Although Rivenbark (2009), found that 
orthorectification was not critical when implementing the shadow technique by-hand, 
FASTEC needs each pixel to be as accurately mapped to the earth as possible, and 
orthorectification was performed. Using fully orthorectified imagery may have 
significantly improved the alignment between files and improved the results (especially 
in complex urban settings). 


21.6 


22.7 



17.4 


16.8 


Figure 15. Nearest neighbor method. In this example, radiance is assigned to a 1 m 
resolution DSM pixel (solid box) using the radiances of four surrounding 2.4m 

satellite pixels (hatched boxes). 
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Due to resolution differences between the DSM and satellite file, weighted 
averaging was used to assign radiances to the DSM pixels (see Figure 15). To begin, the 
latitude (cp D (r)) and longitude (k D (r)) of the DSM pixel at row and column (r,c) are 
determined: 

(p D (r) = cp DuL - (r - l)dcp D ( 13 ) 

1d(c)^ Dul +( c -P d O 4 ) 

where the coordinates of the top left corner ((P Dul , ^ Dul ) anc * change in latitude and 

longitude (dcp D ,d)i D ) of the DSM are read from the DSM flat file header. Next, the 

latitude and longitude are converted into the matching row and column of the satellite file 
(rs, cs) using the following equations: 


r _ (<P SuL -<p D (r)) 11 

dtps 

c _ (^d( c )~^s ul ) | 1 

dl s 


(15) 

(16) 


where the coordinates of the top left comer ((p SlJ[ , ^ S | r ) anc * change in latitude and 


longitude (dcp s ,dA s ) of the satellite file are read from the satellite flat file header. The 


satellite row and column number are real numbers, although row and column indexing 
uses integers. The decimal portion of these numbers is used to detennine the weight of 
each of the four satellite pixels that surround the DSM pixel. Following Franke and 
Nielson (1980), the weights are assigned using inverse weighting based on that decimal 
portion using variations of the following equations (top left equations shown): 



where wtiL is the weight of the top left pixel, R is a constant representing the maximum 
possible fractional distance between the four satellite pixels. Similar equations are used 
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to calculate the weights of the top right (TR), bottom left (BL), and bottom right (BR) 
pixels. The weighted mean of the four satellite pixels relative radiance values (DN T l, 
DN T r, DN B l, DN B r) is calculated using: 

j-jN _ (wt tlDNtl + wt TR DN TR + wt BL DN BL +wt BR DN BR ) 
wt TL +wt TR +wt BL +wt BR 

This value was converted to absolute spectral radiance (in RU) using Equation (8). 

Each satellite image is shipped with a collection of files, including the satellite file 
and a metadata file. These files contain critical infonnation about the location of the 
image, the channel calibration infonnation, and the sun and satellite geometry. A 
procedure was written to read the key data from these metadata files. The IKONOS 
metadata files are fonnatted differently than the files for WorldView-1 and QuickBird. 
They also lacked the necessary calibration infonnation. To conect these issues, a 
program was written to reformat the IKONOS metadata file into a suitable fonn and to 
provide channel-specific calibration information that was obtained from GeoEye 
Corporation. 


3. Output Data 

FASTEC generates two output files to separate general infonnation from the 
AOD calculations and simplify statistical analysis. The first file, identified by a 
“geninfo” tag, contains general information about that particular program run. This 
information includes the satellite name, image date and time, the satellite resolution, 
number of channels, the sun and satellite azimuth and elevation, and the absolute 
calibration factor, effective bandwidth, solar irradiance and 5r for each satellite channel. 
The “geninfo” file also includes the name and resolution of the DSM file, the range of the 
DSM over which the program was run, key shadow/sunlit adjustment variables, the 
number of shadows identified, the DSM offset correction, the number of shadows where 
AOD calculations were performed, and a line stating “Program ended successfully” if the 
program completed without errors. 
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The second output file, identified with “AOD” contains the AOD calculation 
data for each shadow where a calculation was performed. This information includes a 
shadow identification number, a shadow pixel row and column number for the DSM and 
satellite file, the latitude and longitude of that pixel, and pixel counts for the shadow and 
sunlit regions. It also includes the channel number, average shadow and sunlit radiances, 
and the computed TOD and AOD value for each satellite channel. 

If a sub-region is less than or equal to 100 by 100 pixels, the user may create 
additional output files by setting a variable in the configuration file. These additional text 
files store critical array data used in the automated program. They are identified by 
“DSM” and a number as described in Table 2. If sky dome blockage is not tested, file 
number 7 contains zeroes. File numbers 8-11 are updated if the satellite file contains 
those channels. For example, if the satellite file contains one channel, only file 8 would 
be updated while a four channel file would update files 8-11. If additional satellite 
channels were available (not applicable to this study), additional files (numbered beyond 
12) would be generated to store the radiance of those channels. 



Contents of file 

1 

Raw DSM elevation data 

2 

Building height 

3 

Distance from shadow generator with satellite geometry correction 

4 

Same as 3 with pixels that are too close shadow generator invalid 

5 

Same as 3 with shadow/non-shadow edge pixels invalid 

6 

Coherent shadow numbers incorporating invalid pixels files 4 and 5 

7 

Dome block % for each valid shadow/non-shadow pixel 

8-11 

Absolute spectral radiance tor each channel and pixel ~ 


Table 2. List of contents of DSM array output files. 


4. Validation Data 

Aerosol Robotic Network (AERONET) data served as ground truth in this 
research. AERONET is a worldwide collection of robotic sun photometers that measure 
AOD and are accurate to within 0.02 (Holben et al. 1998). Three levels of AERONET 
data are accessible on-line via the Goddard Space Flight Center Web site. Level 1.0 is 
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unscreened data, level 1.5 is cloud-screened data, and level 2.0 is quality assured and 
cloud-screened. Except where otherwise noted, level 2.0 data was used for validation. 

To correct for wavelength differences between AERONET and satellite sensors, 
coefficients for the third order polynomial best-fit curves approximating the spectral 
distribution were determined using Microsoft Excel. The coefficients were determined 
using calculations based on time-averaged AOD values for each wavelength. Although 
Excel has a built-in function for generating polynomial best fit curves, an error in the 
Excel program (detennined by comparing the built-in plot of the curve with a manual plot 
of the curve using the Excel-identified coefficients) produces the wrong coefficients for 
polynomial curves. Therefore, the coefficients were determined using the Excel LINEST 
function (Eng-tips, cited May 2010). These curves, along with the channel-specific 
spectral response curves for each satellite, were used to calculate sensor-specific ground 
truth AOD using Equation (9). IKONOS spectral response data is available on the 
GeoEye Web site while QuickBird and WorldView-1 spectral response curves are 
available from DigitalGlobe upon request. A third order polynomial best fit curve 
estimated AOD values to within 0.003 of truth (see Figure 16). 



300 400 500 600 700 800 900 1000 1100 

Wavelength (nm) 

Figure 16. Third order polynomial best fit curve with raw AERONET data points. 

Four test locations were chosen based on the availability of overlapping 
AERONET stations and Buckeye data. The four stations were Tucson (32N, 110W), MD 
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Science Center (39N, 76W), University of Houston (29N, 95W), and GSFC (38N, 76W). 
Due to the lack of high AOD cases in the available satellite images for the four test 
locations, an additional DSM file was created by hand near the Solar Village (24N, 46E) 
AERONET station (see Figure 17). The heights of the structures were estimated by using 
the length of a shadow and the elevation angle of the sun using the Solar Village image 
from April 12, 2006 at 0754Z. The automated program was designed to calculate AOD 
from the two satellite dish shadows. The building was included to ensure the automated 
method did not include those pixels as part of the sunlit region. The satellite dishes are 
movable objects, and orientation changes did occur between images. Orientation 
differences are noted when observed in testing. 



Figure 17. Satellite image of area used for hand-generated DSM (image captured from 
Google maps). The labeled objects were the only structures built into the DSM. 

D. COMPUTATIONAL TECHNIQUES AND TESTING 

FASTEC was written in FORTRAN with OpenMP for parallel processing. Using 
parallel processing improved the computation time from 18.2 minutes to 4.2 minutes for 
the entire GFSC DSM file using two Intel Xeon 3.2 GHz quad core CPUs running 
LINUX (see Figure 18). The improvements in calculation speed achieved through 
parallel programming made it possible to quickly test key variables across ranges of 
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possible values. Each major step in FASTEC is performed using a different subroutine to 
allow for easy modification of the overall process. The main program performs a few 
basic tasks and calls applicable subroutines as needed. 
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Figure 18. Plot of computation time vs. number of processors. Test case found more 
than 50,000 shadows in 25 km 2 DSM file with AOD calculations on 3015 

shadows. 


The design of the software allowed for pixel-level examination of the automated 
approach in single shadow scenes. Initial testing used this capability to study the impact 
of different shadow and sunlit target selection criteria and led to baseline values for AOD 
computations with every available image and DSM file to test the overall accuracy of 
FASTEC. A separate variable combination was used for the hand-generated Solar 
Village DSM since the shadow generators were comparatively small. For example, the 
satellite dishes in Solar Village were approximately 13 meters tall while nine shadow 
generators in GSFC were at least 35 meters tall. 
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IV. RESULTS 


FASTEC was tested to ensure that it properly identified building shadows and 
calculated AOD. Initial testing examined the primary target selection criteria. The 
results of these tests were used to assign configuration values to perform the next tests 
comparing FASTEC AOD results with ground truth AERONET data. In the case of 
multispectral images, different channels were also compared with each other. Finally, the 
horizontal resolution of the DSM was reduced to examine the coarsest, usable resolution. 

A. TARGET SELECTION CRITERIA TESTS 

The impact of key target selection criteria was analyzed using three satellite 
images; a QuickBird image from October 26, 2006 at 1600Z over the GSFC DSM, a 
WorldView-1 image from January 16, 2009 at 1600Z over the GSFC DSM, and a 
QuickBird image from April 12, 2006 at 0754Z over the hand-generated Solar Village 
DSM. First, pixels along the shadow edge were examined to identify partial shadow 
pixels. After establishing an edge pixel rejection value for the GSFC building, the 
minimum shadow generator height was chosen to prevent AOD calculations on tree 
shadows. This was followed by examining the effect that distance from the shadow 
generator, sky dome blockage, and extreme radiance rejection had on computed values of 
AOD. Except where otherwise noted, the configuration file variables were assigned the 
values listed in Appendix B (Solar Village) and C (GSFC). 

1. Rejection by Minimum Shadow Generator Height (GSFC) 

The minimum shadow generator height was detennined for the GSFC scene after 
the initial single shadow edge pixel test for the GSFC building (see Chapter IV, Section 
A.2.b) detennined the edge pixel rejection variable should be set to 2. For this test, 
FASTEC was run using the October 26, 2006 GSFC image with different values for the 
minimum shadow generator height. Whenever the elevation difference between a 
shadow pixel and the shadow generator was below this threshold, the shadow pixel was 
marked invalid. Next, each shadow in the AOD output file was located in the satellite 
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image by hand and identified as being caused by either trees or buildings. Figure 19 
displays the results of this test with minimum shadow generator height values between 25 
and 40. Lower values were not analyzed by hand due to the large number of shadows 
identified and the fact that 88.2% of the shadows were tree shadows when the minimum 
shadow generator height was set to 25. With settings of 25 or above, the number of 
building shadows only varied by two buildings. With settings of 35 or 40, one tree 
shadow remained. 


■ ■ ■ ■ 

25 30 35 40 

Minimum shadow generator height 

Figure 19. Number of tree/building shadows vs. minimum shadow generator height from 

October 26, 2006 GSFC QuickBird image. 

The variation in AOD values was also examined as the minimum shadow 
generator height was changed (see Figure 20). In this graph, the 50th percentile of AOD 
was plotted for different settings of minimum shadow generator height. The error bars 
represent the 25th and 75th percentile of AOD values, the square marks the 50th 
percentile, the plus symbols mark the maximum and minimum AOD values, and the 
dotted line is the ground truth AERONET AOD value. In cases where the minimum 
shadow generator height is less than 30, the 50th percentile value and range is larger than 
cases where it is 30 or larger. This is due to the large number of tree shadows in the low 
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minimum shadow generator height test cases because tree shadows may allow some 
sunlight to pass through the canopy, increasing the shadow radiance and the computed 
AOD. The maximum value of AOD (not plotted) in test cases below 30 falls between 
0.96 and 2.754. In cases where the number of building shadows exceeds the number of 
tree shadows, the range of AOD values between the 25th and 75th percentile shrinks 
from approximately 0.10 to less than 0.02. In each of these high threshold cases, the 
maximum is due to the one tree shadow remaining in the sample. This tree shadow can 
be eliminated by increasing the minimum pixel threshold to 70. The minimum pixel 
threshold sets the required number of pixels in the shadow and sunlit region for AOD 
calculations to be perfonned. Removing this tree shadow from the sample shrinks the 
AOD range between the maximum and minimum to approximately 0.02, similar to the 
range between the 25th and 75th percentile. Due to these findings, tests involving the 
entire GSFC DSM file used a minimum shadow generator height of 35 and single shadow 
tests on the GSFC building used a minimum shadow generator height of 15 to capture the 
approximately 20 meter tall GSFC building. However, minimum pixel count remained 
15 and any tree shadows were eliminated from the results manually. 
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Figure 20. Minimum shadow generator height vs. 50th percentile of AOD. Computed 
for all shadow generators from the October 26, 2006 GSFC image. Error bars are 
the 25th and 75th percentile, plus symbols are the minimum and maximum AOD 
(maximum exceeds 0.35 in some cases), and the dotted line is ground truth. 
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2 . 


Rejection of Edge Pixels 


Hypothesis I states that it is important to avoid partial shadow pixels because they 
increase shadow radiance and decrease sunlit radiance. The edge pixel rejection variable 
eliminates pixel along the shadow/sunlit edge. Figure 21 shows how the shadow region 
for the north satellite dish (left) and the GSFC building (right) change when the edge 
pixel rejection variable is set to two. In this figure, yellow represents the shadow 
generator, dark grey represents full shadow pixels and light grey represents shadow 
pixels that were eliminated due to their proximity to the shadow edge. The white pixels 
between the shadow and shadow generators are shadow pixels that were eliminated due 
to excessive elevation difference or because they were blocked from view of the satellite. 
For example, in the GSFC scene, pixels on the northwest side of the shadow region were 
eliminated because they were blocked from view of the satellite while some of the pixels 
to the southeast of the shadow region were eliminated due to elevation difference caused 
by a one story annex off that wing of the building. 



Figure 21. Example of edge rejection algorithm on Solar Village north satellite dish (for 
April 12, 2006 image) and GSFC building (for January 16, 2009 image). Dark 
grey pixels are full shadow pixels, light grey are rejected due to proximity to 
shadow edge, and yellow approximate the shadow generator. 


a. North Satellite Dish, Solar Village 

Figure 22 displays the average pixel radiance based on their proximity to 


the shadow edge for the north satellite dish in the April 12, 2006 Solar Village image. 
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The zero line marks the approximate position of the shadow edge so that every pixel falls 
on one side of the shadow edge or the other. The average radiances of the pixels closest 
to the shadow edge fall between more interior pixels. Additionally, the standard 
deviation of pixels that are fully in shadow is low compared to the standard deviation of 
the edge shadow pixels. This is partially due to the different percentages of these pixels 
in direct sunlight but position error between pixels in the DSM and satellite image also 
contributes. Regardless, including these edge pixels increases the range of radiances in 
both the shadow and sunlit region, which should increase the uncertainty of the computed 
AOD from one shadow to the next. For this case, setting the edge pixel rejection variable 
to one eliminates partial shadow pixels without reducing the number of usable pixels 
below the arbitrarily set minimum threshold of 15. Setting the interior rejection variable 
to two offers the lowest uncertainty of the shadow region but it decreases the number of 
usable shadow pixels 62%. 
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Figure 22. Average pixel radiance versus proximity to shadow edge for the north satellite 
dish with the April 12, 2006 Solar Village image. The x-axis is the distance from 

the shadow edge (in pixels) with zero marking the shadow edge, error bars are 
one standard deviation, and plus symbols are the max and min radiance. 
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When calculating optical depth for this shadow, the value decreases as the 
edge pixel rejection variable increases (see Table 3). This decrease is most significant as 
the edge pixel rejection variable changes from 0 to 1. In this example, the computed 
AOD value is close to ground truth when the edge pixel rejection variable is set to 0 
because the partial shadow pixels decrease the radiance difference by over 11 RU when 
compared to the setting of 1. As before, in theory these pixels have some direct sunlight 
reaching them (partial shadow) and the apparent improvement in accuracy is due to 
misapplying the shadow method. For the remainder of tests involving the Solar Village 
scene, the edge pixel rejection variable was set to 1 to minimize the uncertainty of the 
shadow pixel radiance while capturing the most usable shadows across all available 
imagery. The extreme pixel rejection variables (Chapter IV.A.5.b) are designed, in part, 
to compensate for the small number of partial shadow pixels that remain in the sample. 


Edge Pixel 
Rejection 
Setting 

Shadow 

Pixel 

Count 

Non- 
Shadow 
Pixel Count 

Shadow 

Radiance 

Non-shadow 

radiance 

TOD 

AOD 

Error 

0 

168 

1380 

80.98 

152.09 

0.303 

0.242 

-0.001 

1 

98 

1084 

75.44 

158.49 

0.255 

0.194 

-0.049 

2 

37 

775 

70.37 

160.89 

0.225 

0.164 

-0.079 


Table 3. AOD with varying edge pixel rejection setting for the north satellite dish 

in the April 12, 2006 Solar Village image. Ground truth AOD was 0.243. 


b. Goddard Space Flight Center Building 

When the automated method was applied on the GSFC building pictured 
in Figure 21, errors in image geo-referencing and solar azimuth and elevation cause 
minor alignment differences between the DSM and satellite file. While offset correction 
improves translation error between the images, it does not account for rotation or 
distortion errors. In addition, the presence of trees to the southwest of the building 
shadow impacts pixel radiance for more pixels in the sunlit region. All these factors 
make this scene more complex than Solar Village. 

Similar to the Solar Village scene, the average shadow radiance does not 
stabilize until moving three pixels in from the shadow edge (see Figure 23). 
Additionally, sunlit pixel radiance is comparatively low when within two pixels of the 
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shadow edge. Beyond two pixels from the shadow edge the average radiance, standard 
deviation, and difference between the maximum and minimum radiance stabilize in the 
shadow region. Variation in the sunlit region is quite high due to the significant 
contribution of the surface reflectance dependent direct radiance signal to the total signal. 
Based on this examination, it appears the edge pixel rejection variable should be set to 2. 



-11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 


Proximity to shadow edge (pixels) 

Figure 23. Average pixel radiance vs. proximity to shadow edge for the GSFC building 
on January 16, 2009. The x-axis is the distance from the shadow edge (in pixels) 
with zero marking the shadow edge, error bars are one standard deviation, and 
plus symbols are the max and min radiance. 

When computing AOD while varying the edge pixel rejection variable for 
the same scene, fluctuations are quite small (see Table 4). Most notable is the change in 
sunlit radiance when transitioning from a setting of one to two. This comparatively large 
change is due the relatively small sunlit pixel count and the tree line to the southwest. 
The tree line shadow causes the sunlit region to shrink on several sides as the edge pixel 
rejection variable is increased, eliminating a large number of low radiance pixels. This 
supports setting the edge pixel rejection variable to two for the GSFC scene. 
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Edge Pixel 
Rejection 
Setting 

Shadow 

Pixel 

Count 

Non- 
Shadow 
Pixel Count 

Shadow 

Radiance 

Non-shadow 

radiance 

TOD 

AOD 

Error 

0 

924 

357 

16.45 

34.43 

0.124 

0.054 

0.018 

1 

848 

256 

16.16 

35.75 

0.116 

0.046 

0.010 

2 

752 

180 

16.07 

36.77 

0.112 

0.042 

0.006 

3 

668 

137 

15.85 

36.75 

0.109 

0.039 

0.003 

4 

571 

101 

15.2 

36.81 

0.102 

0.032 

-0.004 

5 

460 

69 

14.8 

36.83 

0.098 

0.028 

-0.008 

6 

351 

40 

14.69 

36.57 

0.097 

0.027 

-0.009 

7 

250 

17 

14.6 

34.99 

0.100 

0.030 

-0.006 


Table 4. AOD with varying edge pixel rejection setting for the GSFC building 

using January 16, 2009 image. Ground truth was 0.036. 


Next, AOD values for the whole GSFC DSM file were computed for 
October 26, 2006 using only building shadows with the settings described in Chapter IV, 
Section A.l. The October 26, 2006 image was chosen over the January 16, 2009 image 
because it overlapped the entire GSFC DSM region, providing nine usable shadows 
instead of the four usable shadows in the January 16 image. For these computations, 
offset correction was enabled, the minimum shadow generator height was set to 35, the 
minimum pixel count was set at 15, and tree shadows were manually eliminated from the 
results. Figure 24 displays the results as the edge pixel rejection variable is changed from 
zero to nine. As the edge pixel rejection variable increased from 0, the 50th percentile 
value of AOD decreased. Setting the edge pixel rejection variable to two, as suggested 
by the previous single shadow test, selects all 9 building shadows while keeping the 
range of computed AOD values small. Settings beyond 5 eliminated some buildings 
from the sample. 
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Figure 24. Aerosol optical depth vs. edge pixel rejection setting for October 26, 2006 

GSFC image. The black box is the 50th percentile, the error bars are the 25th and 
75th percentile, the plus symbols are the max and min, and the dashed line is 

ground truth. 

3. Rejection Due to Distance from Shadow Generator 

Hypothesis II states that it is important to choose shadow pixels far from the 
shadow generator to calculate optical depth. Higher elevation objects that surround a 
point block portions of the sky, reducing diffuse radiation components of the total 
radiance signal received at the satellite sensor. Figure 25 shows output from FASTEC 
when a percentage of pixels near the generator are eliminated from the total shadow 
region for both the north satellite dish in Solar Village and the GSFC building. The light 
grey region identifies pixels that are too close the shadow generator. As the percentage 
increases, more of the total shadow region is kept. Testing of the importance of this 
parameter was done using these scenes. Additionally, AOD was computed for the GSFC 
DSM file as the building shadow percentage threshold was varied with 10% increments. 
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Figure 25. Example of target pixel rejection based on proximity to the shadow generator. 
Yellow indicates building pixels, dark grey indicates shadow pixels, and light 
grey indicates eliminated shadow pixels. Created using north satellite dish in 
Solar Village on April 12, 2006 (top panels) and the GSFC building on January 

16, 2009 (bottom panels). 

a. North Satellite Dish, Solar Village 

Figure 26 depicts the spectral radiance based on distance from the shadow 
generator for the April 12, 2006 Solar Village image. For this test case, the edge pixel 
rejection variable was set to 1 to eliminate most of the partial shadow pixels from the 
sample. In general, spectral radiance increases as the distance from the shadow generator 
increases. The large radiance of pixels identified as being within 1.3 meters of the 
satellite dish is heavily influenced by partial shadow pixels in the southeast comer of the 
shadow. The correlation between the distance from shadow generator and spectral 
radiance is about 0.45 when including every shadow pixel, while eliminating the closest 
pixels from the sample increases correlation to 0.77. If the target selection variable that 
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rejects pixels based on their physical distance from the shadow generator were set to 2 or 
rejection by percentage were set to 80%, those partial shadow pixels would be marked 
invalid and eliminated from the calculation. 



Distance from shadow generator (meters) 


Figure 26. Spectral radiance vs. distance from the shadow generator (binned) for the 
Solar Village north satellite dish on April 12, 2006. For each bin, the box marks 
the 50th percentile, error bars are the 25th and 75th percentile, and plus signs are 

the maximum and minimum. 


Table 5 displays the results for AOD calculations on the north satellite 
dish when varying the proximity to building by percentage threshold for the shadow 
region. Generally, AOD increases as more pixels close to the shadow generator are 
removed. The slight AOD decrease between 70% and 80% is due to higher radiance 
from partial shadow pixels near the satellite dish, which are included with the 80% 
threshold. The average shadow radiance increases from 73.04 RU to 78.06 RU as the 
percentage of building shadow is decreased from 80% to 30%. At 80% or greater, the 
entire shadow is used in the computation. Overall, AOD varies by 0.025 across the range 
of thresholds. Fifty percent of the change occurs from 30 to 40% when the pixel count 
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nearly doubles and the average shadow radiance decreases from 78.06 to 75.43 RU. In 


every case, the computed AOD is at least 0.04 lower than ground truth. 


% of shadow 
used 

Shadow 
pixel count 

Non-shadow 
pixel count 

Shadow 

Radiance 

Non-shadow 

Radiance 

AOD 

Error 

100 

98 

1084 

73.04 

159.7 

0.179 

0.067 

90 

98 

1084 

73.04 

159.7 

0.179 

0.067 

80 

98 

1084 

73.04 

159.36 

0.180 

0.066 

70 

86 

1084 

72.89 

159.36 

0.179 

0.067 | 

60 

73 

1084 

73.14 

159.7 

0.180 

0.066 

50 

55 

1063 

74.42 

159.48 

0.187 

0.059 

40 

38 

1021 

75.43 

160.09 

0.191 

0.055 

30 

21 

933 

78.06 

160.47 

0.204 

0.042 


Table 5. Results of varying proximity to building threshold percentage for north 

satellite dish from April 12, 2006 Solar Village image. Ground truth is 0.246. 


b. Goddard Space Flight Center Building 

Figure 27 shows the radiance for pixels binned by distance from the 
shadow generator for the GSFC building on January 16, 2009. In general, the 50th 
percentile value of the radiance increases slightly with distance from the building. The 
high radiance pixels near the shadow generator appear to be partial shadow pixels. Part 
of the variation in spectral radiance across the shadow is due to variations in surface 
reflectance from trees, cement, and grass surfaces in the scene. Other possible reasons 
for the variations across the shadow include rotation or distortion error between the DSM 
and satellite file and reflected direct sunlight off buildings near the shadow. The 
correlation between distance from the shadow generator and spectral radiance for all 
shadow pixels is 0.20, with it increasing to 0.31 when pixels less than 20 meters from the 
building are removed from the sample. This is much lower than the 0.45 and 0.77 
correlations for the Solar Village scene. 
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Distance from shadow generator (meters) 


Figure 27. Spectral radiance vs. distance from the shadow generator (binned) for GSFC 
building on January 16, 2009. For each bin, the box marks the 50th percentile, 
error bars are the 25th and 75th percentile, and plus signs are the maximum and 

minimum. 


Figure 28 is a plot of pixel radiances in the shadow region from the same 
scene. The sun is located to the SSE, casting the shadow to the NW. Subjectively, 
radiance appears to increase with distance from the shadow generator. Additionally, 
radiance is highest along the shadow edge, particularly along the northern edge of the 
shadow. This high radiance region to the north may be due to reflection off the wing of 
the building to the north or west. In both this scene and the Solar Village scene (not 
pictured), there is a region in the center of the shadow that appears to have the lowest 
radiance values. To further examine this relationship, another shadow in the GSFC DSM 
(pictured in Figure 25) was examined. The shadow cast by this tall, isolated building was 
much longer than wide. In this case, the correlation between distance from the shadow 
generator and spectral radiance was only 0.09 (results not pictured). 
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Figure 28. Plot of spectral radiance by pixel for the GSFC building from January 16, 
2009. The GSFC building surrounds all sides except the southwest. The sun is 

located to the southeast. 

Table 6 displays computed AOD values for the GSFC building as the 
building shadow percentage threshold changes from 100 to 20% in 10% increments. For 
this test, offset correction was enabled and the edge pixel rejection threshold was set to 
two. Generally, AOD increased as the threshold decreased and limited the shadow region 
to pixels further from the building. The overall variation in AOD is lower than the Solar 
Village case partly because of the larger number of shadow pixels associated with this 
building (667) and the large width of the shadow generator (about 23 pixels) relative to 
the length (about 15 pixels). Otherwise, the GSFC building and Solar Village scene show 
increasing AOD due to increased shadow radiance as more of the pixels near the building 
are rejected. 
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% of shadow 
used 

Shadow 
pixel count 

Non-shadow 
pixel count 

Shadow 

Radiance 

Non-shadow 

Radiance 

AOD 

Error 

100 

751 

179 

15.65 

37.92 

0.034 

0.002 

90 

751 

179 

15.65 

37.92 

0.034 

0.002 

80 

733 

179 

15.64 

37.92 

0.034 

0.002 

70 

626 

179 

15.74 

37.92 

0.035 

0.001 i 

60 

501 

167 

15.9 

37.82 

0.037 

-0.001 I 

50 

388 

143 

16.11 

37.52 

0.040 

-0.004 

40 

276 

101 

16.41 

37.62 

0.043 

-0.007 

30 

168 

71 

16.77 

38.84 

0.044 

-0.008 

20 

68 

28 

17.28 

44.03 

0.035 

0.001 t 


Table 6. Results of varying proximity to building threshold percentage for GSFC 

building wing from January 16, 2009 GSFC image. Ground truth is 0.036. 


Figure 29 displays the range of AOD values for the entire GSFC DSM 
using the October 26, 2006 image. For this test, offset correction was enabled and the 
building shadow percentage was varied from 10% to 100% with a 10% increment. For 
thresholds above 50%, a lone tree shadow with high AOD (greater than 0.25) was 
eliminated manually. For thresholds 50% and below, tree shadows were not present. 
Otherwise, this graph shows minor variation in both the 50th percentile and range of 
AOD values for all shadows in the scene. In every case, the result shows a low bias of 
approximately 0.05 and the spread across all AOD values varied between 0.026 and 
0.055. These results do not show the percent of building shadow threshold as particularly 
important to the shadow method although the low optical depth of the scene may 
contribute to the lack of change in AOD. For scenes other than Solar Village, the 
building shadow percentage threshold will be set to 50% to increase the average radiance 
slightly. To keep the Solar Village shadow pixel counts above the minimum pixel 
threshold and reduce partial shadow pixels near the shadow generator, the physical 
distance from the shadow generator threshold will be set to two. 
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Figure 29. AOD vs. percent of building shadow used for the October 26, 2006 GSFC 
image. The black box is the 50th percentile, the error bars are the 25th and 75th 
percentile, the plus symbols are the max and min, and the dashed line is ground 

truth. 


4. Rejection Due to Sky Dome Blockage 

Hypothesis III states that sky dome blockage is important when choosing valid 
shadow and sunlit pixels to calculate AOD. This is due to blocked diffuse sky radiance 
(negative contribution to pixel radiance) and direct reflection off taller objects (positive 
contribution to pixel radiance). The sky dome blockage threshold was varied on the 
January 16, 2009 image over the GSFC building to verily the sky dome blockage 
procedure worked properly. Figure 30 displays the change to the shadow region with the 
sky dome blockage threshold set to 30%. In this example, light grey signifies invalid 
shadow pixels, which exceeded the sky dome blockage threshold, dark grey identifies 
valid shadow pixels, and yellow approximates the GSFC building. The invalid region is 
due to higher sky dome blockage from the building wing to the northeast of the shadow 
and to a lesser extent the wing to the northwest. 
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Figure 30. Example of sky dome blockage for the January 16, 2009 image of the GSFC 
building. Yellow is building pixels, dark grey is valid shadow pixels, and light 
grey is invalid shadow pixels from 30% or greater sky dome blockage. 


a. North Satellite Dish, Solar Village 

The north satellite dish in Solar Village provided a simple scene to test the 
significance of sky dome blockage as it relates to spectral radiance. The presence of the 
small building to the west of the dish impacts sky dome blockage and provides a way to 
test the accuracy of the calculation. The maximum vertical blockage of the 4 m tall 
building is about 42 degrees for the closest shadow pixels. The building affects up to 13 
of the 36 azimuths used in the computation with the azimuth interval set to 10 degrees. 
Therefore, the maximum contribution of this building to sky dome blockage within the 
shadow is about 0.119 (see Figure 31). By subtracting the sky dome blockage values 
with the building removed from the blockage with the building included, the maximum 
difference was 0.119 and decreased as distance from the building increased. Considering 
the typical sky blockage of shadow pixels within this shadow is about 0.30, the maximum 
contribution of the building to sky dome blockage is 40%. The correlation between sky 
dome blockage and distance from the shadow generator is -0.79 with the building and - 
0.87 without the building. Correlation is not -1.0 because of the presence of the satellite 
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dish to the south, variation in the north satellite dish pixel elevations, and possible 
incomplete sampling angles in the sky dome blockage calculation. 



Figure 31. Sky dome blockage test using April 12, 2006 Solar Village image. Grey 
region is shadow pixels identified by the automated method. Numbers mark 
geometrically relevant information for the calculation on the highlighted pixel. 

Figure 32 is a graph of the spectral radiance versus sky dome blockage, in 
2.5% bins, for shadow and sunlit pixels for the north satellite dish. Partial shadow pixels, 
because of the edge pixel rejection variable being one, explain the large maximum 
radiance values in some of the shadow bins. Regarding shadow pixels, radiance 
decreases with increasing dome blockage until the 40% bin where the radiance steadily 
increases. The increase may be due to diffuse radiance from the 4 m tall building, partial 
shadow pixels along the shadow/shadow generator edge, or reflection of direct sunlight 
off the side of the building into the shadow. When the sky dome blockage threshold for 
the shadow region was varied from 30% to 50%, AOD changed by less than 0.01 (graph 
not pictured). Dropping the threshold below 30% provided less than the 15 valid shadow 
pixels necessary to perform the calculation. 
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The variation in the radiance of the sunlit pixels in Figure 32 appears 
random with no discernable trends as sky dome blockage increases. This is likely due to 
the overwhelming influence of the direct portion of the total signal radiance. The 
relatively low minima between 15 and 25% dome blockage were due to five partial 
shadow pixels along the edge of the building to the west of the satellite dish. AOD varied 
by less than 0.003 as the sunlit pixel sky dome blockage maximum threshold increased 
from 5% to 100%. Based on these results, sky dome blockage does not appear to be very 

important for either the shadow or sunlit region in this Solar Village image. 
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Figure 32. Spectral radiance of sky dome blockage bins for the April 12, 2006 Solar 
Village image. For each bin, the box (shadow) or diamond (sunlit) is the 50th 
percentile, error bars are the 25th and 75th percentile, plus symbols are the 
maximum, and “x” symbols are minimum radiances. 
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b. Goddard Space Flight Center Building 


The Goddard Space Flight Center building scene provided real DSM data 
to analyze the impact of sky dome blockage on shadow and sunlit regions. The 20 meter 
tall building is composed of several wings off a curved connecting building (see Figure 
33). The shadow region is northwest of the southernmost wing with the sunlit region 
west of the shadow. The least-blocked shadow pixels are in the southwest portion of the 
shadow while the most-blocked shadow pixels are along the northeast edge. The least- 
blocked sunlit pixels are on the southeast side while the most-blocked shadow pixels are 
along the north side of the sunlit region. 



Figure 33. Shadow used for sky dome blockage testing for Goddard Space Flight Center 

building, January 16, 2009 image. 

Figure 34 is a graph of the spectral radiance as sky dome blockage, in 4% 
bin intervals, varies for shadow and sunlit pixels for the GSFC building using the January 
16, 2009 image. For the shadow bins, the maximum radiance pixels fall along the outer 
edge of the shadow suggesting they may be partial shadow pixels caused by position 
error, despite offset correction being performed, between the DSM and satellite file. As 
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sky dome blockage percentage increases, radiance decreases slightly, only to increase 
above 32%. For shadow pixels, correlation between sky dome blockage percentage and 
spectral radiance is 0.35. 

As blockage increases for sunlit pixels, the spectral radiance also increases 
but 144 of the 179 sunlit pixels are less than 16% blocked. The lack of blocked pixels 
above 16% helps explain the large changes in the 50th percentile value of spectral 
radiance at those blocked percentages. Correlation between sky dome blockage and 
spectral radiance is 0.41 in the sunlit region. As discussed previously, when examining 
sunlit pixels, the direct reflection portion of the total signal dominates and it is difficult to 
determine the impact of sky dome blockage. If it were possible to remove the direct 
reflectance portion of the signal, a pattern might emerge. Otherwise, radiances in the 
sunlit region are highest near the building. This is likely due to higher surface reflectance 
for those pixels although additional radiance from reflection of direct sunlight off the 

building may contribute. 

60 

56 

52 

48 

— 44 
ZD 

cc 

<D 40 
o 
c 
gj 

t 36 

aj 

b 32 

CD 
Q_ 

w 28 
24 
20 
16 
12 

0 4 8 12 16 20 24 28 32 36 40 44 48 52 

Sky dome blockage (4% bin interval) 

Figure 34. Spectral radiance of sky dome blockage bins for the GSFC building on 
January 16, 2009. For each bin, the box (shadow) or diamond (sunlit) is the 50th 
percentile, error bars are the 25th and 75th percentile, plus symbols are the 
maximum, and minus symbols are the minimum radiances. 
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Table 7 displays the results of FASTEC as the shadow sky dome blockage 
maximum threshold is varied. Although the AOD values are quite small in each case, the 
average shadow radiance and AOD increase as the sky dome blockage threshold 
increases. Part of the increase is due to the higher radiances associated with high sky 
dome blockage pixels although this increase is partially countered by variation in offset 
correction as the threshold changes (see Table 9). The change in offset correction values 
also contributes to the variation in sunlit radiance. If offset correction perfectly aligned 
the DSM and satellite files, sunlit radiance would vary much less than it does in this 
table. 


&naaow sxy 
dome blockage 
threshold 

Shadow 
pixel count 

Non-shadow 
pixel count 

Shadow 

Radiance 

Non-shadow 

Radiance 

AOD 

Error 

0.15 

20 

152 

14.09 

34.31 

0.026 


0.20 

184 

179 

14.42 

34.66 

0.028 

-0.008 

0.25 

324 

180 

14.75 

37.6 

0.025 

-0.011 

0.30 

411 

180 

14.81 

37.6 

0.026 

-0.010 

0.35 

509 

180 

15.03 

38.39 

0.026 

-0.009 

0.40 

614 

179 

15.06 

37.72 

0.028 

-0.008 

0.45 

708 

179 

15.36 

37.92 

0.031 

-0.005 

0.50 

750 

179 

15.65 

37.92 

0.034 

-0.002 

0.55 

751 

179 

15.65 

37.92 

0.034 

-0.002 


Table 7. Results of varying shadow sky dome blockage maximum threshold for 

January 16, 2009 image of GSFC building wing. Ground truth is 0.036. 


Table 8 shows the results from FASTEC as the sunlit sky dome blockage 
maximum threshold is varied for the GSFC building shadow using the January 16, 2009 
image. Generally, sunlit radiance increases as the sunlit sky dome blockage threshold is 
increased as expected from Figure 34. Since shadow radiance is fixed, AOD decreases 
but the overall change is quite small. This is due to the majority of sunlit pixels having 
sky dome blockage below 15%. 
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Non-snaaow sKy 
dome blockage 
threshold 

Shadow 
pixel count 

Non-shadow 
pixel count 

Shadow 

Radiance 

Non-shadow 

Radiance 

AOD 

Error 

0.05 

751 

22 

15.65 

36.87 

0.037 

0.001 | 

0.10 

751 

96 

15.65 

37.37 

0.035 

-0.001 

0.15 

751 

137 

15.65 

36.66 

0.037 

0.001 

0.20 

751 

156 

15.65 

36.93 

0.036 

0.001 

0.25 

751 

169 

15.65 

37.42 

0.035 

-0.001 

0.30 

751 

173 

15.65 

37.52 

0.035 

-0.001 

0.35 

751 

178 

15.65 

37.83 

0.034 

-0.002 

0.40 

751 

179 

15.65 

37.92 

0.034 

-0.002 | 

0.45 

751 

179 

15.65 

37.92 

0.034 

-0.002 


Table 8. Results of varying the sunlit sky dome blockage maximum threshold for 

January 16, 2009 image of GSFC building wing. Ground truth is 0.036. 


Finally, the sky dome blockage thresholds were varied for the entire GSFC 
DSM file using the October 26, 2006 image. Figure 35 is a plot of the results of running 
FASTEC while varying the sunlit sky dome blockage maximum threshold. This plot 
shows little variation in AOD calculations as the threshold changes. This is due to the 
fact that almost all the pixels in the sunlit regions have sky dome blockage less than 0.05 
as evidenced by the average sunlit pixel count increasing from 563 to 677 as the 
threshold increases from 0.05 to 0.30. At thresholds beyond 0.30, sunlit pixel count does 
not change. 
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Figure 35. AOD vs. sunlit pixel sky dome blockage maximum thresholds for the GSFC 

DSM file from January 16, 2009. For each threshold, the box is the 50th 
percentile, the error bars are the 25th and 75th percentile, and plus symbols mark 

the maximum and minimum AOD. 

Figure 36 is a plot of the results of running FASTEC while varying the 
shadow sky dome blockage maximum threshold. The minimum AOD (-0.11) is below 
the y-axis for the 0.05 threshold. Setting shadow sky dome blockage thresholds to more 
than 0.60 does not eliminate any pixels. As opposed to changing sunlit dome blockage 
maximum threshold, despite small AOD values, changes in the 50th percentile values are 
apparent. These changes mirror changes in the row and column offset for the different 
thresholds (see Table 9). Where the row and column offset values are similar, the results 
in the figure are similar. For example, when the threshold is 30% or higher, offset 
correction adjustments are within two pixels of the baseline and the difference between 
the maximum and minimum AOD is less than 0.002. Based on the inconclusive results 
of these tests, sky dome blockage thresholds will be disabled for later testing. 
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Figure 36. AOD vs. shadow pixel sky dome blockage maximum thresholds for the GSFC 
DSM file from January 16, 2009. Boxes are the 50th percentile, the error bars are 
the 25th and 75th percentile, and + symbols mark the max and min AOD. 


Shadow sky 
dome blockage 
threshold 

Row 

Offset 

Col 

Offset 

Distance from 
baseline (2,9) 

0.05 

13 

11 

11.18 

0.10 

10 

9 

8.00 

0.15 

10 

9 

8.00 

0.20 

9 

9 

7.00 

0.25 

6 

8 

4.12 

0.30 

4 

9 

2.00 

0.35 

3 

8 

1.41 

0.40 

3 

8 

1.41 

0.45 

3 

9 

1.00 

0.50 

1 

8 

1.41 

0.55 

2 

9 

0.00 


Table 9. Offset correction values for different shadow sky dome blockage 

thresholds for January 16, 2009 image over the GSFC DSM. Baseline offset 
correction determined by disabling the sky dome blockage maximum threshold. 
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5. 


Rejection of Extreme Radiances 


a. North Satellite Dish, Solar Village 

Hypothesis IV states that eliminating the extreme radiances values from 
the sample improves AOD calculation. Figure 37 shows the sorted spectral radiances of 
shadow and sunlit pixels for the north satellite dish in Solar Village from the April 12, 
2006 image using FASTEC. The first approximately 100 sunlit pixels increase in 
radiance steeply due to partial shadow pixels along the northwest side of the small 
building west of the satellite dish and partial shadow pixels from both satellite dish 
shadows. Including these pixels in the sunlit region average causes a slight decrease in 
the average radiance for the region although the large number of sunlit pixels lessens the 
influence of these darker pixels. When the 105 darkest (lowest 10%) pixels are removed 
from the sample, the average sunlit radiance increases by 1.4 RU, however, when the 
darkest and brightest extremes are removed equally across the sunlit region, the change in 
radiance is quite small (see Table 11). 

For shadow pixels, the overall slope of the increasing radiances is steeper 
and the six brightest pixels have noticeably higher radiances than the rest of the sample. 
The overall steeper slope of the shadow region when compared with the sunlit region 
appears to be partly due to the smaller number of pixels in the shadow region (98 vs. 
1059 pixels). The significantly brighter pixels are partial shadow pixels along the 
shadow edges. Removing these six pixels causes average radiance to decrease by almost 
1.2 RU. Removing these extreme radiance pixels from both regions results in a 2.6 RU 
increase in the radiance difference and computed AOD to be .0083 lower. Although this 
lowering of AOD pushes the results further from ground truth, it eliminates invalid sunlit 
pixels and properly applies the shadow method. 

When the darkest and brightest extremes are removed equally from the 
shadow region, the largest change in shadow radiance occurs when 10% of the extremes 
are removed (see Table 10) and radiance difference increases from 86.8 RU to 87.6 RU, 
decreasing AOD from 0.192 to 0.171. The radiance difference change is less than 0.2 
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RU as additional percentages of the extremes are eliminated. For sunlit pixels, radiance 
difference changes by 0.3 or less across all combinations (see Table 11). 



+ Non-shadow 
- Shadow 


Sorted pixels, Solar Village 04122006 

Figure 37. Spectral radiance of usable pixels for shadow and sunlit regions for Solar 
Village north satellite dish and the April 12, 2006 satellite image. 


EH 

EM 

Min pixel 
position 

Max pixel 
position 

Min pixel 
radiance 

Max pixel 
radiance 

Average 

shadow 

radiance 

Radiance 

difference 

1 o 

100 

1 

98 

66.3 

96.8 

73.0 

86.8 1 

10 

90 

9 

88 

67.9 

77.9 

72.3 

87.6 

20 

80 

19 

78 

68.9 

75.7 

72.2 

87.6 

30 

70 

29 

68 

70.3 

74.0 

72.2 

87.7 

40 

60 

39 

58 

71.1 

73.4 

72.1 

87.7 

50 

50 

49 

49 

72.1 

72.1 

72.1 

87.8 


Table 10. Results of varying extreme radiance thresholds for Solar Village north 
satellite dish shadow region and the April 12, 2006 satellite image. 
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EH 

EH 

Min pixel 
position 

Max pixel 
position 

Min pixel 
radiance 

Max pixel 
radiance 

Average non¬ 

shadow 

radiance 

Radiance 

difference 

0 

100 

1 

1059 

115.3 

179.3 

159.9 

86.8 

10 

90 

105 

953 

152.8 

168.8 

160.1 

87.0 

20 

80 

211 

847 

155.0 

165.6 

159.9 

86.9 

30 

70 

317 

741 

156.8 

163.0 

159.8 

86.8 

40 

60 

423 

635 

158.4 

161.2 

159.8 

86.7 

50 

50 

529 

529 

159.7 

159.7 

159.7 

86.7 | 


Table 11. Results of varying extreme radiance thresholds for Solar Village north 
satellite dish sunlit region and the April 12, 2006 satellite image. 


b. Goddard Space Flight Center Building 

Figure 38 is a plot of the sorted spectral radiance of usable shadow and 
sunlit pixels for the GSFC building using the January 16, 2009 image. In this case, the 
shadow region is much larger than the sunlit region and the subsequent slopes are 
reversed with the smaller sunlit region (179 pixels) being steeper than the shadow region 
(751 pixels). Otherwise, the shape of the shadow region curve is similar to the shadow 
region curve from Solar Village with the highest radiance pixels (32 of 748) being 
significantly brighter than the rest of the dataset. Interestingly, the sunlit region curve has 
a different shape than Solar Village with the brightest pixels being much brighter than the 
rest of the sample due to variations in surface reflectance from two distinct surface types; 
some pavement with a primarily grassy surface. Otherwise, the impact on the radiance 
difference as larger percentages of the extremes from the total sample are removed is 
very similar to the Solar Village case where removing 10% of the darkest and brightest 
pixels has the largest affect on the average radiances (results not pictured). 

Although this suggests that 10% is the ideal threshold of the extremes, the 
small variation as more of the extremes are rejected leaves open the possibility that larger 
portions of the extremes can be trimmed without significantly changing the resulting 
AOD values. The benefit to trimming larger percentages of the extreme radiances should 
be most evident when calculating AOD over heterogeneous surfaces or when overcoming 
issues with offset correction or partial shadow pixel rejection thresholds. The results 
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from the GSFC and Solar Village tests both indicate that any threshold over 10% should 
provide representative samples of the shadow/sunlit regions. 



+ Non-shadow 
- Shadow 


Figure 38. Spectral radiance of usable pixels for shadow and sunlit regions for GSFC 

building and the January 16, 2009 image. 


Figure 39 displays the results from running FASTEC on the GSFC DSM 
for the October 26, 2006 image with different thresholds of extreme radiance rejections. 
In each case, both extremes were trimmed an equal percentage for both the shadow and 
sunlit regions. This graph highlights the fact that computed AOD values are very 
consistent for each of the nine usable shadows in this image, regardless of the threshold 
setting. In this scene, 3 of the 9 AOD values are consistently negative regardless of the 
threshold. These negative values occur when molecular Rayleigh optical depth is 
subtracted from total optical depth. Otherwise, computed AOD values for each building 
vary by less than 0.02 as the extreme rejection values change. These facts suggest the 
threshold has little impact on the results, however, the true AOD for this scene is 0.059, 
well below the recommended threshold of 0.10 (Vincent 2006). In addition, the 
minimum building height setting happens to select buildings whose shadows fall on 
largely homogeneous surfaces, reducing the impact of these variations of the extreme 
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rejection threshold. Based on the findings from this section, the extreme rejection 
thresholds will be set to 25% for the whole DSM tile tests in the next section. This 
ensures the middle 50% of radiances in the shadow/sunlit regions are used to compute 
optical depth and should minimize variations in AOD, making the results more precise. 
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Figure 39. Aerosol optical depth using different extreme radiance rejection percentages 
on the GSFC DSM file and the October 26, 2006 image. For each threshold, the 
box is the 50th percentile, the error bars are the 25th and 75th percentile, plus 
symbols mark the maximum and minimum, and the dotted line marks ground 

truth AOD. 

B. WHOLE SCENE AOD CALCULATION 
1. Solar Village 

The first full scene investigation of FASTEC was conducted using the manually 
created Solar Village elevation file and 14 satellite images (5 multispectral, 9 
panchromatic) with corresponding ground truth AERONET data. This scene consisted of 
simple elevation data with high surface reflectance in a desert region. Figure 40 is the 
manually created 64 cm horizontal resolution elevation file covering an approximately 
3600 square meter area about 200 meters northeast of the AERONET station (located at 
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N 24° 54’ 25” N, E 46° 23’ 49”). The surface elevation was set to 1 m, the rectangular 
building elevation was set to 5 m, and the satellite elevations were set to about 14 meters. 
Appendix D lists the configuration settings with minimum building height set to six to 
ensure AOD calculations were only performed on the two satellite dishes, not the small 
rectangular building. The satellite dishes were movable and their orientation was 
occasionally different than pictured. 



Figure 40. Manually created DSM for Solar Village. The (1) north and (2) south satellite 
dishes were estimated to be 13 meters tall and the building was estimated to be 4 

meters tall. 

a. Multispectral Images 

Figure 41 displays the red, green, and blue (RGB) channels from the five 
multispectral QuickBird images that were evaluated using FASTEC. The NIR channel 
was also evaluated (not pictured). Ground truth best fit curves were determined with 
level 2.0 AERONET daily average AOD whenever it changed by less than 0.02. In the 
case of the April 30, 2006 image, level 2.0 AERONET data stopped at 0707 and the 
image was collected at 0755. Level 1.0 data showed AOD increasing rapidly during the 
satellite collection time. Surface visibility at Riyadh International Airport, located 
approximately 15 miles southeast of the test location, decreased from 8 km to 4 km 
between 0800 and 1000 (the 0900 observation is missing). With a westerly wind of 8-10 
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knots during that period, it would take between one and two hours for dust to advect over 
Riyadh from the region where the AERONET site is located. Based on Koschmieder’s 
equation (Horvath 1971), these AOD values are consistent with observed surface 
visibility (see Equations (20) and (21)): 


Vis = 


3.912 

5 „ 


a„ = 


H 


( 20 ) 

( 21 ) 


where 5 a is AOD, a e is the extinction coefficient, and H is the height of the planetary 


boundary layer. For this day, the mean of the two surrounding level 1.0 AERONET 
samples was used instead of the daily average. In the October 19, 2005 and the 
November 30, 2007 images, satellite dish orientations were noticeably different from the 
elevation file. 



Figure 41. Multi-spectral imagery for Solar Village testing on (A) October 19, 2005, (B) 
April 12, 2006, (C) April 30, 2006, (D) May 26, 2007, and (E) November 30, 
2007 with 2.4 meter resolution QuickBird data. 


Basic data for each image test is listed in Table 12. Relatively large 
shadows (with a pixel count of more than 100) were identified in the October and 
November images while the April and May images had relatively small shadows. Offset 
correction shifted the DSM less than 20 meters in any direction for each image. In the 
May 26 image the mean sun azimuth and elevation information in the metadata file that 
accompanied the image was incorrect (150.2 and 44.9 respectively) and was adjusted by 
hand based on historical data from the Naval observatory. 
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Date 

Time 

Sun 

Azimuth 

Sun 

Elevation 

Row 

Offset 

Column 

Offset 

Satellite 

Dish 

Shadow 
Pixel Ct 

Non-shadow 
Pixel Ct 

10/19/2005 

0752 

160.2 

53.4 

-17 

31 

1 

113 

643 

2 

101 

658 

4/12/2006 

0754 

135.4 

68.4 

3 

-2 

1 

27 

551 

2 

33 

588 

4/30/2006 

0755 

124.3 

73.5 

-8 

7 

1 

13 

476 

2 

11 

576 

5/26/2007 

0801 

105.2 

77.6 

-4 

-13 

1 

6 

382 

2 

9 

553 

11/30/2007 

0754 

164.9 

42.2 

-2 

29 

1 

107 

469 

2 

101 

560 


Table 12. Date, time, sun azimuth, offset correction, and pixel counts for the satellite 
dish shadows in the Solar Village multispectral imagery. 


In general, AOD values using FASTEC were quite accurate with bias less 
than 0.05 and standard deviation less than 0.08 (see Figures 42-45 and Table 13). Using 
Equations (20) and (21), this standard deviation equates to a horizontal visibility variation 
of less than 1.5 km when AOD is 0.5 or greater and the planetary boundary layer height 
is 1 km. The largest error was found in the two April images. On April 30, 2006, level 
1.0 data showed rapidly increasing AOD to over 1.0 around the time the image was 
collected. It is interesting that FASTEC showed an unusually large decrease in AOD 
between channels when comparing April 30, 2006 with the other days. This may be due 
to the size of the aerosols present or due to problems with the shadow method when 
optical depth is high. Otherwise, the 10 data points showed a slight high bias for the blue 
channel and a slight low bias for the others with each channel error standard deviation 
less than 0.08. 
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Green channel AOD (FASTEC) Blue channel AOD (FASTEC) 



■ 10192005 
♦ 04122006 
a 04302006 
X 05262007 
+ 11302007 


Figure 42. Blue channel AOD results for Solar Village satellite dishes. 



■ 10192005 
♦ 04122006 
a 04302006 
X 05262007 
+ 11302007 


Figure 43. Green channel AOD results for Solar Village satellite dishes. 


72 
















Near-Infrared channel AOD (FASTEC) Red channel AOD (FASTEC) 


1.2 


0.8 


0.6 


0.4 



■ 10192005 
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Figure 44. Red channel AOD results for Solar Village satellite dishes. 



■ 10192005 
♦ 04122006 
a 04302006 
X 05262007 
+ 11302007 


Figure 45. Near-Infrared channel AOD results for Solar Village satellite dishes. 
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Channel 

Bias 

Standard 

Deviation 

Blue 

0.043 

0.078 

Green 

-0.006 

0.078 

Red 

-0.008 

0.079 

NIR 

-0.003 

0.078 


Table 13. Multispectral channel-specific error bias and error standard deviation 

Solar Village shadows. 


b. Panchromatic Images 

Nine panchromatic images over Solar Village were evaluated with 
FASTEC (see Figure 46). The satellite dishes are pointing in different directions in 
several of the images, most notably (F) and (I). This should affect the accuracy of the 
shadow found using the hand-generated Solar Village DSM with the dish orientation 
pointing almost straight up. Ground truth best fit curves were based on the daily average 
AERONET value for each image except for the April 30, 2006 image. As discussed in 
the previous section, level 1.0 AERONET data showed AOD increasing rapidly from 
about 0.7 to 1.0 around the time of the satellite image. For every other image, the AOD 
changed slightly with standard deviations less than 0.02 for the whole day. 
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Figure 46. Panchromatic imagery for Solar Village testing on (A) October 19, 2005, (B) 
April 12, 2006, (C) April 30, 2006, (D) May 26, 2007, (E) November 30, 2007, 
(F) June 26, 2008, (G) October 11, 2008, (H) October 14, 2008, (I) December 14, 

2008. 

As in the multispectral images, FASTEC computed AOD for both satellite 
dishes in each image. Basic data on each of the scenes is provided in Table 14. Four 
images, collected between April and June, had relatively small shadows with fewer than 
33 pixels while the others identified between 80 and 200 pixels. The horizontal 
resolution of WorldView-1 data (50 cm) was slightly higher than DSM resolution (64 
cm) while the QuickBird imagery had the same resolution as the DSM. As in the 
multispectral case, the May 26, 2007 sun azimuth and elevation was incorrectly identified 
in the metadata file and was manually corrected to 105.2° and 77.6° respectively. 
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Date 

Sat 

Time 

Sun 

Azimuth 

Sun 

Elevation 

Row 

Offset 

Column 

Offset 

Satellite 

Dish 

Shadow 
Pixel Ct 

Non-shadow 
Pixel Ct 

10/19/2005 

QB 

0752 

160.2 

53.4 

-14 

31 

1 

113 

643 

2 

101 

658 

4/12/2006 

QB 

0754 

135.4 

68.4 

4 

-1 

1 

27 

551 

2 

33 

588 

4/30/2006 

QB 

0755 

124.6 

73.5 

-38 

14 

1 

13 

477 

2 

11 

577 

5/26/2007 

QB 

0801 

105.2 

77.6 

-13 

-2 

1 

6 

382 

2 

9 

553 

11/30/2007 

QB 

0754 

164.9 

42.2 

4 

30 

1 

108 

469 

2 

101 

560 

6/26/2008 

WV1 

0736 

90.3 

71.5 

2 

10 

1 

18 

446 

2 

20 

475 

10/11/2008 

WV1 

0743 

154.6 

55.1 

-7 

17 

1 

94 

628 

2 

82 

635 

10/14/2008 

QB 

0750 

158.6 

54.9 

-25 

-4 

1 

93 

614 

2 

83 

646 

12/14/2008 

WV1 

0746 

161.3 

39.7 

4 

22 

1 

199 

505 

2 

167 

488 


Table 14. Date, satellite, time, sun azimuth and elevation, offset correction, and 

pixel counts for the satellite dish shadows in the Solar Village panchromatic 

imagery. 


Figure 46 is a plot comparing the AOD values found using the satellite 
images with AERONET. The results showed a low bias of 0.136 and an error standard 
deviation of 0.071 for all 18 AOD calculations. If the suspect April 30, 2006 case is 
removed from the sample, the low bias becomes 0.115 and the error standard deviation 
lowers to 0.042. Figure 48 displays the bias corrected results. Although the sample is 
rather small and the bias and accuracy are found using the one scene, this figure shows 
good agreement between ground truth and FASTEC. 
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■ 04122006 
a 04302006 
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x 06262006 
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05262007 


Figure 47. Panchromatic channel AOD results for Solar Village satellite dishes. 



♦ 10192005 
■ 04122006 
a 04302006 

11302007 
x 06262006 

• 10112008 
+10142008 
- 12142008 

05262007 


Figure 48. Panchromatic channel bias-adjusted AOD results for Solar Village satellite 

dishes. 
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2 . 


Goddard Space Flight Center 


The first full scene investigation of FASTEC using a BuckEye DSM was 
conducted with the approximately 5 km x 5 km Goddard Space Flight Center (GSFC) 
file. The region provides a mixture of urban and rural terrain with large buildings, 
asphalt, and forested areas (see Figure 49). Appendix E lists the configuration values 
where the minimum building height and minimum pixel count were set up to eliminate 
tree shadows. Using these settings prevented AOD calculations on all but nine buildings 
from the DSM, including the Goddard Space Flight Center building where the 
AERONET station was located and where single shadow testing in Chapter IV, Section B 
was performed (see upper right hand corner of Figure 49). 



Figure 49. BuckEye DSM tile over the Goddard Space Flight Center with zoomed in 
view of three buildings identified as shadow generators by FASTEC. Green dots 
are shadow regions with sun to SE. Red dot marks location of AERONET site. 
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a. 


Multispectral Images 


Two multispectral images that covered the GSFC region and coincided 
with AERONET data were evaluated with FASTEC. Figure 50 shows the RGB images 
for the zoomed in region of Figure 49. The fourth channel, at the NIR wavelength, was 
also evaluated (not pictured). Ground truth best fit curves were based on the daily 
average AERONET value for each image because AOD changed less than 0.02. Portions 
of the AERONET data was missing from the cloud screened data on February 8, 2009. It 
is easy to identify the building shadows in the October 26, 2006 image on the left while 
the shadows in the February 8, 2009 image on the right are covered with thick cumulus 
clouds. 



Figure 50. Examples of multispectral imagery for GSFC BuckEye file testing from 
October 26, 2006 (left) and February 8, 2009 (right). Images are approximately 
equivalent to the zoomed in region in Figure 49. 


On February 8, 2009, every usable target shadow is obscured by cloud. 
For 4 of the 6 target shadows, FASTEC found shadow radiance exceeded sunlit radiance 
and AOD was not calculated. The other target shadows calculated unrealistic AOD 
values, between 0.43 and 1.13 with ground truth of 0.07 for the red channel for instance. 
For this reason, the February 8, 2009 image was not included in the results. 

On October 26, 2006 at 1621Z, the sun azimuth and elevation were 170.7° 
and 38.2° respectively. Offset correction shifted the DSM file eight meters east to align it 
with the satellite pixels. For the eight usable shadows, the average shadow region was 
273 pixels while the average sunlit region was 296 pixels. The 50th percentile value of 
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AOD for the blue channel was slightly above ground truth while the 50th percentile value 
was below ground truth for the other channels (see Figure 51). For the red channel, the 
calculated AOD was less than zero for each of the eight shadows and the 50th percentile 
was 0.091 below ground truth (see Table 15). For each channel, the standard deviation 
was 0.078 or better. For the blue channel, the bias was noticeably above ground truth due 
to one shadow, which found AOD to be 0.295. Removing that shadow from the dataset 
brought the bias to within 0.003 of ground truth and reduced the standard deviation to 
0.051. 



♦ Blue 

■ Green 
▲ Red 

■ NIR 


Figure 51. 50th percentile value of AOD vs. AERONET value for GSFC multispectral 
channels from October 26, 2006 QuickBird image. Error bars on the y-axis 
reflect the reported accuracy of AERONET data (0.02) and error bars on x-axis 
are one standard deviation of the values reported on that day and scene. n=8. 
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Channel 

Bias 

Standard 

Deviation 

Blue 

0.027 

0.078 

Green 

-0.055 

0.023 

Red 

-0.091 

0.039 

NIR 

-0.018 

0.020 


Table 15. Multispectral channel-specific error bias and standard deviation for 

October 26, 2006 GSFC shadows. 

Interesting results appear when plotting all other channels against the blue 
channel (see Figure 52). It is expected that AOD would decrease with increasing 
wavelength, however, FASTEC AOD values decreased from blue to green to NIR to red. 
The two exceptions were the cases where blue channel AOD was lowest. In the case 
where the green and blue channel had the lowest AOD values, the sunlit region included 
water. In the other case, green was lower than NIR but higher than red and the sunlit 
region was composed of pavement, grass, and cars. 

The maximum and minimum blue channel AOD cases are also noticeably 
different than the other six cases. These extremes were primarily due to different surface 
types dominating the sunlit region. In the maximum blue channel AOD case, the sunlit 
region surface consisted primarily of grass with two sidewalks and a dirt area present. 
This resulted in sunlit radiance being significantly lower than the six pavement cases 
(36.5 RU vs. 44-50 RU). In the minimum blue channel AOD case, the sunlit region 
surface was primarily composed of sidewalk and water with some grass as well. This 
resulted in sunlit radiance being significantly higher than the six pavement cases (58 RU 
vs. 44-50 RU). In both the blue channel extremes, the green and NIR AOD was 
consistent with the other shadows. For the maximum blue channel AOD, the red channel 
AOD was much lower than the other shadows because its sunlit radiance was 
significantly lower (19 RU vs. 25-30 RU). Overall, shadow radiance varied by about 4 
RU for each channel except NIR where it varied by 8 RU. 
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Figure 52. Plot of blue channel AOD vs. the other channels for the October 26, 2006 

GSFC image. 


b. Panchromatic Images 

Figure 53 shows the four panchromatic images that were evaluated with 
FASTEC. Surface obscuring clouds were not apparent for any of the valid target 
shadows in the four satellite images although thin cirrus may have been present. Ground 
truth was determined using best fit curves based on the daily average AERONET value 
for each image because AOD changed by less than 0.02 for each wavelength. As 
explained in the multispectral section, cloud coverage over every target shadow in the 
February 8, 2009 led to unusable results and the image was discarded. 


82 










Figure 53. Examples of panchromatic imagery for GSFC BuckEye file testing from (A) 
October 26, 2006, (B) March 5, 2008, (C) May 25, 2008, and (D) January 16, 
2009. Images are approximately equivalent to the zoomed in region in Figure 49. 


In each image, FASTEC computed AOD for anywhere from three to eight 
shadows, depending on the overlap between the DSM file and the satellite image. Basic 
data is listed in Table 16 where the top image was from QuickBird, the bottom three from 
WorldView-1. Generally, shadow size was largest for the January 16, 2009 image where 
the sun elevation was lowest. Offset correction shifted the DSM between about 3 and 14 
meters. The sun azimuth and elevation provided in the metadata files appeared accurate 


when subjectively compared with the satellite images. 


Date 

Time 

Sun 

Azimuth 

Sun 

Elevation 

Row 

Offset 

Column 

Offset 

Shadow 

Count 

Avg. Pix 

Shadow 

el Count 

Non-shadow 

10/26/2006 

1621 

170.7 

38.2 

2 

9 

8 

273 

296 

3/5/2008 

1543 

148.1 

40.4 

-2 

3 

7 

215 

308 

5/25/2008 

1606 

140.8 

68.4 

-1 

-13 

4 

79 

332 

1/16/2009 

1600 

159.6 

27.7 

-7 

5 

3 

560 

365 


Table 16. Date, time, sun azimuth and elevation, offset correction, shadow count, 

and average pixel counts for the satellite dish shadows in the Solar Village 

panchromatic imagery. 
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Figure 54 displays the results for the four test cases. In each case, AOD 
values were less than 0.1, with FASTEC finding the lowest values for the May 25, 2008 
image. The March 5 and January 16 images showed the 50th percentile AOD value 
slightly above ground truth while the other two images were below ground truth. The 
standard deviation of the AOD values of the May 25 image was slightly less than 0.03 
while the other images were less than 0.02. Overall, the four test cases suggested a low 
bias of approximately 0.04. In the May 25, 2008 case, the negative AOD values occurred 
after removing the Rayleigh contribution. 



■ 20061026 

♦ 20080305 
■ 20080525 

♦ 20090116 


Figure 54. 50th percentile value of AOD vs. AERONET value for GSFC panchromatic 
images. Error bars on the y-axis reflect the reported accuracy of AERONET data 
(0.02) and error bars on x-axis are one standard deviation of the values reported 

for that day and scene. 


3. Other Location Tests 


a. Tucson 

Tucson was chosen to test FASTEC using BuckEye DSM data in a desert 


region. Figure 55 shows the DSM file used for this location with the zoomed in region 

84 












showing the area surrounding the AERONET station. The configuration values were the 
same as those used for GSFC testing (see Appendix E). The large bright region in the 
southwest corner of the DSM is caused by mountainous terrain. 



Figure 55. DSM file for Tucson testing. The inset is the area around the AERONET 
station, the bright area in the southwest corner is caused by mountains. 

Only one available Tucson image corresponded with AERONET data. 
The image for that day, November 1, 2005, was comprised of the panchromatic and blue 
channels from the IKONOS satellite (see Figure 56). The four meter resolution blue 
channel is much fuzzier than the one meter resolution panchromatic channel. Ground 
truth was determined using best fit curves of the daily average AERONET AOD because 
it varied by less than 0.02 for each wavelength. 
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Figure 56. November 1, 2005 1815Z imagery for Tucson. Image (A) is the four-meter 
resolution blue channel and (B) is the one-meter resolution panchromatic channel. 

The area in each image roughly matches the zoomed in area in Figure 55. 

FASTEC computed AOD for both channels. The sun azimuth and 
elevation were 163° and 41.6° respectively. Offset correction adjusted the DSM 13 rows 
and 24 columns for the blue channel and 14 rows and 24 columns for the panchromatic 
channel. Calculations were perfonned on 14 shadows in the blue channel with an 
average shadow region size of 257.3 pixels and average sunlit region size of 256.4 pixels 
while 3 shadows in the panchromatic channel had an average size of 199.3 shadow pixels 
and 214.7 sunlit pixels. The disparity in the number of targets between the two channels 
was because the areal extent of the panchromatic channel was much smaller than the blue 
channel. 

The blue channel AOD values had a significant high bias of over 0.28 (see 
Figure 57). The standard deviation of AOD values for the blue channel was 0.055, 
suggesting the method found values in good agreement with each other despite being 
very high. It appears the high bias was caused by the poor resolution collecting data over 
a larger area for each pixel in the satellite image. This combined with the weighting 
method used to assign radiance values to the DSM pixels smoothed the pixel radiances, 
decreasing the radiance difference between shadow/sunlit regions, and causing AOD to 
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be too high. The complex city scene may amplify this problem and IKONOS 

multispectral imagery may work properly in scenes with more homogeneous surfaces. 
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Blue channel error bin (computed - truth) 

Figure 57. Histogram of blue channel AOD computation error for November 1, 2005 

image. Ground truth AOD was 0.046. 

The panchromatic channel expected AOD had a slight low bias of 0.029 
relative to ground truth of 0.026, with each computed AOD value being slightly less than 
zero (-0.0045, -0.0021, -0.0008). The standard deviation of AOD values for the 
panchromatic channel was 0.001, showing strong agreement. The findings for this 
channel support the belief that the comparatively poor resolution of the IKONOS 
multispectral sensor may be insufficient for this scene. 

b. Baltimore and Houston 

Baltimore and Houston were included to test FASTEC in very complex 
settings with numerous shadow generators located in close proximity to each other (see 
Figure 58 and 59 respectively). Of the three Houston images tested, none had coinciding 
ground truth AERONET data while 14 of the 15 Baltimore images coincided with 
AERONET data. Due to problems with FASTEC in these complex settings, comparisons 
with ground truth were not perfonned. 
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Figure 58. Baltimore DSM. The zoomed in region shows good shadow generators near 

the Maryland Science Center AERONET station. 



Figure 59. Houston DSM. The zoomed in region shows good shadow generators near 

the University of Houston AERONET station. 
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The Houston DSM included the downtown area with a complicated 
shadow patterns (see Figure 60). Although FASTEC appeared to correctly identify 
shadow regions, the offset correction algorithm failed to correctly align the DSM and 
satellite file. This led to widely varying, inaccurate AOD values. For example, the 
January 25, 2004 image found 24 targets for AOD calculation. 3 of the 24 had shadow 
radiance greater than sunlit radiance and AOD was not calculated. The other 21 targets 
had an average AOD in the blue channel of 0.338 and a standard deviation of 0.449. 
Although ground truth was not available, the comparatively large standard deviation 
suggested the results were invalid. When the locations of the shadows in the output file 
were compared with the satellite image, it was obvious offset correction failed to 
correctly align the satellite and DSM. Attempts to align them by hand were not 
successful either. The other Houston images suffered from poor offset correction as well. 



Figure 60. Multispectral QuickBird image of Houston downtown area from January 25, 
2004 at 1647Z. Note the complexity of shadows cast by the tall buildings. 


89 


The Baltimore DSM included the downtown area and Chesapeake Bay 
(see Figure 61). The large body of water further complicated the scene and contributed to 
the failure of offset correction and AOD values to be wrong. It is theorized that the low 
radiance over water offered a large area where the difference between shadow and sunlit 
radiance was quite small, contributing to the failure of the offset correction algorithm to 
properly adjust the DSM. For the January 1, 2004 panchromatic image, the average 
AOD was 0.305 and the standard deviation was 0.345 suggesting the results were invalid. 
When the location of shadows in the output file were compared with the actual shadows 
in the satellite image, it was apparent that offset correction again failed. In fact, offset 
correction failed to correctly align the DSM and satellite file in every Baltimore and 
Houston image although it worked properly in images from Solar Village, GSFC, and 
Tucson. Several attempts to manually correct for the offset failed to adequately align the 
images. As a result, the results from complex urban scenes were discarded for this study. 



Figure 61. Panchromatic QuickBird image of Baltimore from January 1, 2004 at 1540Z. 
Note the complexity of shadows cast by the tall buildings and the optically dark 

Chesapeake Bay. 
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4. 


Summary of Whole Scene AOD Calculations 


FASTEC showed promise when used in the Solar Village and GSFC scenes. In 
both scenes, small bias pointed to problems with accuracy but FASTEC calculated AOD 
values that were consistent with each other, as evidenced by standard deviations less than 
0.08. Results from the GSFC scenes suggest that bias correction may work best when the 
surface type is known. The test of the IKONOS image in Tucson suggested the lower 
resolution multispectral sensor may be problematic, at least in urban settings, but the 
panchromatic channel appeared to work properly. The complex scenes in Houston and 
Baltimore highlighted problems with the offset correction procedure and identified the 
need to improve upon that procedure or use higher quality orthorectified satellite imagery 
before testing FASTEC in such complex urban settings. 

C. DSM RESOLUTION ADJUSTMENT 

The final test of FASTEC was designed to determine the minimum usable 
resolution of elevation data by reducing the horizontal resolution of the DSM. The 
horizontal resolution of the DSM was reduced using two methods. Method one sampled 
a subset of the full resolution DSM data while method two averaged the full resolution 
DSM data (see Figure 62). The first method is designed to simulate what the DSM might 
look like if a LIDAR sampled the surface with half the frequency of BuckEye while the 
second method was designed to simulate what the results might look like using stereo 
pairs DSM data. 
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Figure 62. Example of resolution adjustment methods. Box A represents a full resolution 
DSM, box B and C halve the horizontal resolution using method one (sampling) 
and two (smoothing) respectively. Gray shading is by elevation. 


1. Resolution Reduction for Solar Village Scene 

Figure 63 shows examples of reduced resolution DSM images using the sampling 
and smoothing method with resolution reduced to 1.28 and 1.92 meters. The most 
notable difference between the two methods is the edges of the objects are less sharp 
when the smoothing method is applied. The most notable difference between the 
resolution changes with either method is the larger area of the full DSM that is visible in 
the 30 pixel by 30 pixel region. Otherwise, resolution reduction using sampling makes 
the edges appear sharper while smoothing makes the edges less sharp. The configuration 
values for this test are listed in Appendix D. Resolution reductions beyond 1.92 meters 
lowered the pixel count below the user-defined minimum of 5. 
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Figure 63. Examples of top left 30 x 30 pixel region in Solar Village DSM after 

resolution adjustment methods applied. 


Table 17 displays the results of resolution reduction testing using the sampling 
method with the Solar Village elevation file and the multispectral satellite channels. As 
elevation fde resolution is reduced, AOD values decrease with the overall decrease 
typically 0.05 for the blue channel and about 0.03 for the green, red, and NIR channels. 
These decreases are primarily due to decreasing shadow radiance. For example, blue 
channel radiance values decrease by about 1 RU for each unit of resolution decrease on 
dish 1 and 3 RU for each resolution change on dish 2. These cause an increase the 
radiance difference and decrease in computed AOD values. The fact that the changes are 
consistent for each shadow suggests lower resolution elevation data may still be usable, 
with only the bias needing adjustment. 
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Date 

Dish 

Channel 

AOD w/ Resolution 

0.64 

1.28 

1.92 

Truth 

10192005 

1 

Blue 

0.301 

0.296 

0.256 

0.372 

Green 

0.241 

0.232 

0.215 

0.347 

Red 

0.223 

0.214 

0.195 

0.317 

NIR 

0.220 

0.212 

0.198 

0.297 

2 

Blue 

0.293 

0.268 

X 

0.372 

Green 

0.253 

0.232 

X 

0.347 

Red 

0.254 

0.231 

X 

0.317 

NIR 

0.254 

0.231 

X 

0.297 

11302007 

1 

Blue 

0.263 

0.247 

0.213 

0.168 

Green 

0.202 

0.181 

0.164 

0.149 

Red 

0.171 

0.155 

0.154 

0.124 

NIR 

0.157 

0.140 

0.144 

0.104 

2 

Blue 

0.266 

0.236 

X 

0.168 

Green 

0.217 

0.193 

X 

0.149 

Red 

0.196 

0.165 

X 

0.124 

NIR 

0.181 

0.156 

X 

0.104 


Table 17. Results of Solar Village resolution reduction using sampling approach and 
the multispectral channels. x=insufficient pixels for calculation. 


Table 18 displays the results of resolution reduction testing using the sampling 
method with the Solar Village elevation file and the panchromatic satellite channel. With 
resolution of 1.92 meters, AOD was not calculated on the southern satellite dish (Dish 2) 
and on the north dish on October 11, 2008. Otherwise, AOD varied less than 0.01 for 
every image except the November 30, 2007 image where it varied about 0.023. Variation 
this small suggests FASTEC works with elevation files of up to 1.92 meter resolution and 
possibly beyond since FASTEC stopped calculating AOD before the values changed 
significantly. 

Results from the smoothing method (not pictured) were very similar to the 
sampling method with one exception. On October 14, 2008 with resolution at 1.92 
meters, the Dish 1 AOD value was 0.777 because shadow radiance increased from 
around 55 RU to 121.4 RU with only 3 pixels. Otherwise, as with the multispectral case, 
the AOD value changes appear consistent for each shadow, suggesting using lower 
resolution data may only affect the bias of the results. 
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Date 

Dish 

AOD w/ Resolution 

AERONET 

AOD 

0.64 

1.28 

1.92 

10112008 

1 

0.095 

0.093 

X 

0.267 

2 

0.096 

0.091 

X 

10142008 

1 

0.113 

0.109 

0.108 

0.245 

2 

0.115 

0.107 

X 

10192005 

1 

0.175 

0.167 

0.165 

0.319 

2 

0.183 

0.174 

X 

11302007 

1 

0.088 

0.084 

0.066 

0.124 

2 

0.087 

0.087 

X 

12142008 

1 

0.091 

0.090 

0.088 

0.230 

2 

0.104 

0.097 

X 


Table 18. Results of Solar Village resolution reduction testing using sampling 

approach and the panchromatic channel. x=insufficient pixels for calculation. 


2. Resolution Reduction for GSFC Scene 

Figure 64 shows examples of reduced resolution DSM images using the sampling 
and smoothing method with resolution reduced to 2 and 3 meters. Although more 
difficult to discern than in the Solar Village image in Figure 63, the edges of the objects 
are less sharp when the smoothing method is applied. This is most noticeable with the 
road running north-south and the trees south of the building. The configuration file was 
set to the values listed in Appendix E with the exception of the minimum pixel count 
variable that was reduced to from 100 to 6. This change was made to allow the resolution 
reduction algorithm to work with the expected lower pixel counts. In cases where this 
change caused a tree shadow to have AOD computed, the tree shadow was removed from 
the sample. 
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2 m 3 m 



Figure 64. Examples of top left 201 x 201 pixel region in GSFC DSM after resolution 

adjustment methods applied. 

Table 19 displays the results of resolution reduction testing on the October 26, 
2006 GSFC multispectral image. The 50th percentile AOD values tend to decrease with 
reduced elevation file resolution, although the change is on the order of 0.02 for all 
channels. These decreases are primarily caused by decreasing shadow radiance that 
increases the radiance difference between the regions. As in the Solar Village test, the 
fact that the changes are consistent for each shadow suggests lower resolution elevation 
data may still be usable, with only the bias needing adjustment. This test suggests 4 
meter horizontal resolution elevation files may be adequate for computing AOD, 
although lower resolution may work with very large objects in simple scenes. 
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DSM 

Number 
of calcs 

Bli 

50% 

je AOD 
Range 
(75%-25%) 

Gre 

50% 

sen AOD 
Range 
(75%-25%) 

Re 

50% 

d AOD 
Range 
(75%-25%) 


m 

Full 

8 

0.121 

0.049 

0.031 

0.017 

-0.017 

0.027 

0.020 

0.024 

Sample 2 

6 

0.107 

0.016 

0.025 

0.011 

-0.017 

0.010 

0.017 

0.016 

Sample 3 

5 

0.106 

0.045 

0.021 

0.007 

-0.022 

0.006 

0.012 

0.011 

Sample 4 

5 

0.098 

0.004 

0.017 

0.015 

-0.037 

0.013 

0.003 

0.008 

Sample 5 

2 

0.084 

0.019 

0.007 

0.002 

-0.031 

0.004 

0.000 

0.007 

Smooth 2 

6 

0.110 

0.021 

0.026 

0.010 

-0.020 

0.023 

0.018 

0.010 

Smooth 3 

5 

0.095 

0.027 

0.016 

0.003 

-0.030 

0.004 

0.017 

0.010 

Smooth 4 

5 

0.088 

0.028 

0.014 

0.004 

-0.029 

0.009 

0.012 

0.003 

Smooth 5 

1 

0.095 


0.002 


-0.050 


-0.016 



Truth 0.104 0.083 0.058 0.040 


Table 19. Results of October 26, 2006 GSFC multispectral resolution reduction 

testing. 
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V. CONCLUSIONS 


A. SUMMARY 

FASTEC worked well for Solar Village and the Goddard Space Flight Center 
testing but failed to work properly for other test locations. Shadow selection criteria 
testing suggested it is important to remove partial shadow edge pixels from the sample. 
This testing also suggested it is somewhat important to eliminate pixels that are close to 
the building and eliminate a portion of the extremes from the shadow and sunlit regions. 
The results from sky dome blockage testing were inconclusive. Whole tile testing does 
suggest it is useful to calculate AOD over a region although the test set-up and limited 
data set prevented fully assessing the accuracy (bias correction) or precision (variation) of 
FASTEC. Finally, the DSM resolution test showed that it is possible to use lower 
resolution elevation files and still achieve good results. Further testing may reveal that 
the elevation file resolution must only be capable of discerning shadows to employ 
FASTEC. 

Strengths of FASTEC include the speed with which this method identifies good 
targets, the objective detennination of shadow/sunlit regions, the easily configurable 
nature of the algorithm, and the consistency in AOD values between targets in the same 
scene. When offset correction worked properly, the technique consistently calculated 
AOD values within 0.05 for targets in a scene. Elevation data resolution appears usable 
up to 4 meters depending on the scene being analyzed. This technique can fill the high 
surface reflectance region gap in remotely sensed AOD databases. Militarily, FASTEC 
could enable in-scene correction in intelligence collection and possibly influence electro- 
optical targeting decisions. 

The primary weakness of this study is the failure of the offset correction 

algorithm in complex urban settings such as Houston and Baltimore. When the output of 

FASTEC has high standard deviation, the user should locate several target shadows to 

determine if the cause of the high standard deviation was due to offset correction error, 

the presence of clouds, or some other problem. Additionally, this study suggested the 
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shadow technique has bias, which may vary between channels, scenes, or surface types. 
Although only a few images were used in this study, the automated approach makes 
future, larger studies using different sensors feasible for shadow analysis. These studies 
should refine the bias and uncertainty of the shadow method. 

FASTEC is capable of supporting operational customers or scientists. An 
intelligence analyst might use FASTEC to determine the optical depth of a scene and use 
that value to perform in-scene corrections to the image. A climatologist or numerical 
modeler might establish a global grid of shadow generators to feed a global aerosol 
database while the ability to examine detailed shadow characteristics allows researchers 
to determine different shadow characteristics. 

B. FUTURE RESEARCH 

Future research should focus on two areas: testing FASTEC with large datasets or 
different sensors and improving the technique itself. When it comes to large datasets, test 
locations should include AERONET data or another reliable, accurate means of 
measuring ground truth. Elevation data can come from LIDAR (as was used in this 
study), stereo pair imagery, or other sources capable of providing 4 meter horizontal 
resolution. Using large datasets may enable the user to sort results by surface type, 
building size, or sensor. The ability to define small sub-regions within the DSM should 
simplify this type of study. If the user analyzed each shadow generator in one full scene 
run and created configuration files around the desired target shadows, the output could be 
sorted for straightforward statistical analysis. Similar experiments employing different 
sensors with different wavelength bands or resolution are also possible. For example, 
hyperspectral images provide hundreds of usable wavelength bands and may even allow 
aerosol size characterization. Additionally, it may be possible to use lower resolution 
satellites exploiting mountain or cloud (if cloud top heights are known) shadows. 
Finally, testing FASTEC in complex urban settings with fully orthorectified satellite 
imagery would determine the feasibility of using FASTEC in such scenes. 

There are several possible ways of improving FASTEC to include eliminating 
dark surfaces, applying a cloud mask, identifying isolated buildings, and improving offset 
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correction. Optically dark surfaces such as water or trees pose a problem for FASTEC 
because they could cause offset correction to fail to correctly align the elevation and 
satellite files. If pixels in the elevation file were identified by surface type, these regions 
could be avoided. Similarly, it is necessary to employ the shadow technique in cloud-free 
regions. Adding a cloud mask would eliminate cloud-obscured targets from the 
calculation and improve the technique. This mask should be applied after identifying 
shadow pixels but before performing offset correction. A robust cloud mask should 
account for the shadow cast by the cloud and block off a standard size region around the 
cloud to eliminate the indirect reflective affects of the cloud. Also, Vincent (2006) 
suggested selecting isolated shadows for AOD calculation. It may improve the technique 
to eliminate target shadows when other shadow generators within a user-defined area are 
present. 

Offset correction should improve when working with fully orthorectified imagery. 
In addition, it may improve with the addition of a dark surface and cloud masks but 
additional changes may be required. First, the current method works by trying to align 
the shadow regions in the DSM with low radiance pixels in the satellite file. 
Incorporating a dark surface mask may eliminate this issue by marking these low 
radiance pixels invalid while a cloud mask would prevent high radiance cloud pixels 
from skewing the correction. If these masks fail to correct the problem, an alternative 
approach is to adjust the DSM position to maximize the radiance difference between 
shadow and sunlit pixels. This approach would take both shadow and sunlit radiance into 
account unlike the current approach that just examines shadow radiance. Additionally, it 
may be better to perform offset correction on each individual shadow instead of the entire 
elevation file. This approach would increase computation time but make each AOD 
calculation independent, improving the likelihood that some, or most, of the shadows are 
properly aligned. Regardless, correcting for the offset appears to be the biggest obstacle 
to implementing this algorithm operationally. 
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APPENDIX A: FORTRAN CONFIGURATION VARIABFES 


A. GENERAL USE VARIABLES (CITY_VARS) 


Name 

Description 

Sat_path 

Path to satellite file 

DSM_path 

Path to DSM file 

op_path 

Path to output directory 

Sat name 

Base input satellite filename 

DSM name 

Base input DSM filename. Do not include a “.hdr” file descriptor 

Op name 

Base output filename. Makes “AOD.txt” and “.geninfo.txt” files 

DSM_sr/er, 

DSM_sc/ec 

Limits shadow search to sub region of DSM file. If set outside 

the bounds of the DSM file, the entire DSM file is explored. 


B. SHADOW REGION VARIABLES (SHADOW_VARS) 


Name 

Description 

DSM min bldg hgt 

Minimum height of shadow generator 

DSM_bldg_dist 

Percentage of shadow near building that is rejected 

DSM min bldg dist 

Minimum distance of valid shadow pixel from shadow generator 

rlus/rlls 

Upper/lower bounds of sorted shadow pixel radiances 

max ediff 

Maximum range of elevation differences among shadow pixels 


C. SUNLIT REGION VARIABLES (SUNLIT_VARS) 


Name 

Description 

DSM u rgn 

Size of region around shadow to search for sunlit pixels 
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nonshad int 

Azimuth interval to use when searching for sunlit pixels 

nonshad dint 

Distance interval on each azimuth for finding sunlit pixels 

nonshad ediff 

Max elevation difference of sunlit pixel and shadow 

rlun/rlln 

Upper/lower bounds of sorted sunlit pixel radiances 


D. SHADOW/SUNLIT REGION VARIABLES (SHADNON_VARS) 


Name 

Description 

DSM how int 

Depth of shadow edge to mark invalid due to edge contamination 

minpix 

Minimum count of shadow/sunlit pixels to calculate AOD 


E. MISCELLANEOUS VARIABLES (MISC_VARS) 


Name 

Description 

whichAz 

Which azimuth (min/mean/max) in metadata file to use 

whichEl 

Which elevation (min/mean/max) in metadata file to use 

edgeadj 

Defunct 

off cor 

Set to .true, to perform automated offset correction 

roff, c off 

Used to manually offset DSM row and column 

offshad 

Defunct 

offset 

Size of adjustment region for the DSM file 

offchan 

Offset correction sat channel (set to 1 if value exceeds # of bands) 

chk 

Defunct 

RR 

Weighting factor for nearest neighbor determination 

DSM invalid 

Value to use for program identified invalid pixels 

bigdatadump 

Set to .true, to dump DSM arrays to text files for detailed analysis 
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defunct 


parallel 


F. VARIABLES FOR DOME BLOCKAGE CALCULATION (DOMEJVARS) 


Name 

Description 

error check 

Set to .true, to perform dome blockage calculation on each pixel 

domesize 

Size of region around pixel to calculate dome blockage 

aziinc 

Azimuth increment for dome calculation 

distinc 

Distance increment for dome calculation 

maxdomes 

Maximum allowed dome blockage for shadow pixels 

maxdomen 

Maximum allowed dome blockage for sunlit pixels 


G. AOD CALCULATION VARIABLES 


Name 

Description 

interval 

Scattering angle interval in MAR calculation 

ssa 

Single-scatter albedo for AOD calculation 

g 

Henyey-Greenstein constant 
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APPENDIX B: LIST OF BASELINE USER-DEFINED VARIABLES 
FOR SHADOW CHARACTERIZATION TESTS (SOLAR VILLAGE) 


0.0 

DSM min bldg hgt 


F 

off cor 

1.0 

DSM_bldg_dist 


0 

r off 

0.0 

DSM min bldg dist 


0 

coff 

1.0 

Rlus 


100 

offshad 

0 

Rlls 


10 

offset 

1.5 

max ediff 


1 

offchan 

15 

DSM u rgn 


T 

bigdatadump 

3 

nonshad int 


T 

parallel 

0.99 

nonshad dint 


T 

error check 

1.0 

nonshad ediff 


20 

domesize 

1.0 

rlun 


10 

aziinc 

0.0 

rlln 


2 

distinc 

1 

DSM how int 


1.0 

maxdomes 

15 

minpix 


1.0 

maxdomen 

2 

whichAz 


10 

interval 

3 

whichEl 


0.88 

ssa 

200 

edgeadj 


0.65 

g 
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APPENDIX C: LIST OF BASELINE USER-DEFINED VARIABLES 
FOR SHADOW CHARACTERIZATION TESTS (GSFC BUILDING) 


15.0 

DSM min bldg hgt 


F 

off cor 

1.0 

DSM_bldg_dist 


0 

r off 

0.0 

DSM min bldg dist 


0 

coff 

1.0 

Rlus 


100 

offshad 

0 

Rlls 


10 

offset 

1.5 

max ediff 


1 

offchan 

15 

DSM u rgn 


T 

bigdatadump 

3 

nonshad int 


T 

parallel 

0.99 

nonshad dint 


T 

error check 

1.0 

nonshad ediff 


20 

domesize 

1.0 

rlun 


10 

aziinc 

0.0 

rlln 


2 

distinc 

2 

DSM how int 


1.0 

maxdomes 

15 

minpix 


1.0 

maxdomen 

2 

whichAz 


10 

interval 

3 

whichEl 


0.88 

ssa 

200 

edgeadj 


0.65 

g 
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APPENDIX D: LIST OF USER-DEFINED VARIABLES FOR 
SOLAR VILLAGE AOD CALCULATIONS 


6.0 

DSM min bldg hgt 


T 

off cor 

0.8 

DSM_bldg_dist 


0 

r off 

2.0 

DSM min bldg dist 


0 

coff 

0.75 

Rlus 


100 

offshad 

0.25 

Rlls 


40 

offset 

1.5 

max ediff 


4 

offchan 

15 

DSM u rgn 


F 

bigdatadump 

3 

nonshad int 


T 

parallel 

0.99 

nonshad dint 


F 

error check 

1.0 

nonshad ediff 


20 

domesize 

0.75 

rlun 


10 

aziinc 

0.25 

rlln 


2 

distinc 

1 

DSM how int 


1.0 

maxdomes 

5 

minpix 


1.0 

maxdomen 

2 

whichAz 


10 

interval 

3 

whichEl 


0.88 

ssa 

200 

edgeadj 


0.65 

g 
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APPENDIX E: LIST OF USER-DEFINED VARIABLES FOR 
BUCKEYE FILE AOD CALCULATIONS 


35.0 

DSM min bldg hgt 


T 

off cor 

0.5 

DSM_bldg_dist 


0 

r off 

2.0 

DSM min bldg dist 


0 

coff 

0.75 

Rlus 


100 

offshad 

0.25 

Rlls 


150 

offset 

1.5 

max ediff 


4 

offchan 

15 

DSM u rgn 


T 

bigdatadump 

3 

nonshad int 


T 

parallel 

0.99 

nonshad dint 


F 

error check 

1.0 

nonshad ediff 


20 

domesize 

0.75 

rlun 


10 

aziinc 

0.25 

rlln 


2 

distinc 

2 

DSM how int 


1.0 

maxdomes 

100 

minpix 


1.0 

maxdomen 

2 

whichAz 


10 

interval 

3 

whichEl 


0.88 

ssa 

200 

edgeadj 


0.65 

g 
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